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ABSTRACT 


The purpose of this thesis is to model a Combat System utilizing modern 
methods of nonlinear nonequilibrium statistical mechanics. This initiates development 
of methods which eventually can be used as a decision aid to the commander in force 
planning, battle management, budgeting decisions, doctrinal evaluations, and combat 
analysis. A general method is developed and then applied to a particular battle 
scenario using the combat wargame JANUS. The method fits empirical data to a 
functional form (a Lagrangian) which describes the short time probability distribution 
of a set of order parameters. A maximum likelihood fit is obtained using a simulated 
annealing optimization algorithm. The most likely states of the order parameters and 
the associated risks (variances) of those states ultimately depend on the detailed 
structure of the Lagrangian. A long time probability distribution of the order 


parameters can then be found by using the path integral. 
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I. INTRODUCTION 


Imagine you are the commander of a large force faced with the following 
situation. You have been ordered to defend a key piece of terrain. Intelligence sources 
estimate that an enemy force, approximately three times as large, is moving towards 
your position and is expected to arrive within a couple of hours. You now must make 
a decision on how to allocate your forces on the line of defense in order to repel the 
enemys attack. You have several alternatives. You could leave your forces as is on a 
line defense. But vou know the enemy will only attack a small portion of your front 
and use his overwhelming force to penetrate your position. You could place them in a 
dispersed defense. You ask your operations officer to develop other alternatives. You 
must have them quickly so that the defense plan can be promulgated to the units in a 
timelv fashion. 

The operations officer and his plans/analysis officer, armed with PIACA 
(pronounced PI-CA), Path Integral Algorithm for Combat Analysis, begin to develop 
the alternatives. PIACA 1s a hardware/software package designed to give the 
commander the most likelv results of decisions and the risks associated with that 
decision. By inputing information about their own forces and those of the enemy, and 
information relating to the tvpe of mission, PIACA will give them the most probable 
outcomes of the forces (levels) at the end of some pre-selected time interval. By 
modifving the scenario slightly as to initial force levels and other parameters, they will 
then have a good idea of the best alternatives to present to the commander. The 
commander has now an objective evaluation of his alternatives and is able to make a 
more informed decision. 

There will be occasions when the commander is under a severe time constraint 
and must make a decision based on incomplete information. He now has PIACA 
available as a powerful aid to combat planning and analysis. It is the purpose of this 
thesis to develop a stochastic model of combat and a generalized methodology based 
on that model for eventual use in PIACA. It is additionally shown how the model and 
the methodology can be used for a simple scenario based on data from the combat 
wargame JANUS. PIACA can also be used to evaluate new system hardware,1.e. 


Weapons systems, analyze combat plans, and aid in the analysis of doctrinal changes. 


I] 


Chapter 2 outlines the Lanchester theory of Combat systems. This chapter 1s 
provided as background. Chapter 3 introduces the concept of order parameters and 
discusses their relation with combat systems. Several possible order parameters for the 
combat system are presented. Chapter 4 develops the underlying mathematical theory 
and introduces the reader to the Lagrangian and the associated Path Integral, a 
mathematical physics approach to C? systems developed by Ingber [Ref. 1,2]. Chapter 
5 develops PIACA as a generalized methodology for modeling combat systems. 
Chapter 6 gives several empirical examples. The first example uses a one order- 
parameter model with simulated data generated from a stochastic Lanchester equation 
With constant variance. This will be shown to be equivalent to a quadratic Lagrangian 
With the result that the distribution of the order parameters will be Gaussian with non- 
stationary means and variances. The second example is a two order parameter model 
using simulated data from a different stochastic Lanchester equation. The long time 
conditional distribution will be shown to be non-Gaussian even though the short time 
distribution is Gaussian with respect to the temporal changes of order parameters. 
These examples are provided to show the relationship between the stochastic 
Lanchester representation and the Lagrangian representation. The third example will 
begin with a Lagrangian representation. Data from the combat wargame JANUS is 
used to develop the functional form (Lagrangian). Then an analysis of the short time 
probability distribution of the order parameters using the Lagrangian is given. Chapter 
7 concludes the thesis with a summary of the methodology, its importance and utility, 


and discusses further applications of PIACA and development of the full decision aid. 


Il. AN INTRODUCTION TO LANCHESTER THEORY 


In this chapter, we outline the Lanchester model of warfare, both deterministic 
and stochastic. For a more detailed development, the interested reader should refer to 
Taylor [Ref. 3}. 


A. DETERMINISTIC MODELS 


Lanchester originally formulated his model of combat as a set of differential 


equations, one being, 


X 


dX/dt = -aY where X(tp) = Xp 


Y 


dY/dt = -bX where Y(t)) = Yo (2.1) 


where X and Y are the number of combatants for each side and a,b are called attrition 


rate coefficients. This is Lanchester’s aimed fire model. The other is 


-ax Y 


w 


0 
| 


-bXY , (2.2) 


where the variables are defined above. Equation 2.1 when integrated yields the so- 


called “square-law” 
b(X)* - X?) = a(Y,~ - Y?) (2.3) 


which gives the interpretation that the more initial force level a side has the less his 
casualties. This equation assumes the casualty rate is proportional to the number of 
combatants. It is also referred to as a state equation. Equation 2.2 1s referred to as 
the state equation for area fire, and assumes fire is distributed uniforniy by the 


combatants. When integrated, equation 2.2 yields 


b(X, - X(t)) = a(Yq - Y(t) (2.4) 


i. 


and is called the Lanchester linear law. Although the above equations model simple 
combat quite well, they are limited in scope and have several disadvantages [Ref. 3]: 

1. Constant rate coefficients 

2. No force movement 

3. Various aspects of logistics, C°I, terrain, geographics, etc. are not considered. 
This has led to various modifications which attempt to overcome the above 
shortcomings such as: 

1. using variable rate-coefficients 

2. modeling breakpoint or battle-termination conditions 

3. modeling morale 

4. modeling communications [Ref. 4], etc. 
The basic disadvantage with using deterministic Lanchester equations is that in teat 


combat is a severely stochastic system, which leads us to our next topic. 


B. STOCHASTIC MODELS 

In attempting to define a stochastic model from a deterministic model we must 
recognize that this definition is not unique, i.e., there is a many-to-one mapping of 
deterministic systems into stochastic contexts. For example, one might arbitrarily add 


noise to equation 2.1 in the following manner, 


X 


(X+n) (2.5) 


@e 


X 


I 


f{X)+f'(X)n+ higher order terms (2.6) 


where ny has some distribution and zero mean. Another example might be to add noise 


in an additive way, such as 


x 


-aY + gn (2.7) 


-< 
| 


= -bX + en (2.8) 
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where g is some constant multiplving the standard deviation of the noise. We could 
also formulate a model stochastically by developing a set of Kolmogorov equations. 
Once developed, all models should be able to answer questions concerning the outcome 
of the battle and other factors such as: 

1. What is the probability of win for each side? 

2. How does win probability change with initial force levels? 
3. What is the expected length of battle? 

4. What is the probability distribution of the force levels? 
As 1s evident, there are a number of possibilities of stochastically modeling combat. In 
this thesis we will develop a stochastic model of combat which, as a side benefit, will 
incorporate an underlying physical explanation. It will have several advantages: 

1. to model the stochastic nature of combat 

2. to answer questions such as those above concerning the battle 


to incorporate non-linearities in the model 


{y2 


= 


with the methodology developed, to be able to fit empirical data to the model 
and thus have the potential of forecasting battle outcomes. 


ie 


Ill. ORDER PARAMETERS AND COMBAT 


A. INTRODUCTION 

We will begin this discussion with assumptions and definitions. A battle will be 
defined as a combat engagement between two opposing forces constrained to a small 
geograpluc area. This will be our system that we will attempt to model. The state of 
the system will be defined as a collection of variables which, as a set describe the 
system at any time, t. 

This system will be nonlinear, dynamic and stochastic: nonlinear, since the 
moments of the distribution of the state variables may be described as nonlinear 
functions of the other state variables; dynamic, since the state variables could be 
functions of time; and stochastic, because the variables will be random due to inherent 
noise in the system. This noise will reflect imperfect knowledge of the enemy’s forces, 
weather, equipment failure, etc., and also of the commander’s own forces, and may also 
be a nonlinear function of the state variables. 

When modeling this system we have several alternatives. One alternative would 
be to use generalized stochastic differential equations as our model with the variables 
denoting the microscopic state of the system. This is a very intuitive approach, but 
there are mathematical difficulties in solving such large sets of coupled stochastic 
differential equations, and even more difficulties in interpreting the numerical results. 
However, there may be alternative sets of variables which define the system 
appropriately enough for a study of combat at a middle, i.e. intermediate level of 
aggregation, or “mesoscopic” level. A probability distribution of these new variables 
would allow us to make predictions of the variables at any time, t. We will call these 
new variables order parameters (Ref. 5]. The order parameters of the system will 
contain all the information inherent in the system, relevant to a specific “coarse- 
grained” scale to be studied, and should be kept at a minimum to allow easy 
assimilation by the commander essentially defining the appropriate scale of aggregation 
to be considered. This is one of the problems associated with complex C3] and combat 
systems; i.e. We must be careful not to pass on too much information to overburden 


the command nodes. 
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B. REPRESENTATIONS OF COMBAT 

We could describe this battle in several different, equivalent representations. For 
example, consider the grid shown in Figure 3.1 as a “coarse” geographic representation 
of our battle. The ier represent tanks of force pt in cell i where ft=1,2 and i= 1-9. 
Sinular notation exists for the personnel, BE . The state variables would then 
incorporate all microscopic information such as velocity, position, ammo, etc. Another 
could be a time sequence event description, ie. at each moment in time, a particular 
event occured and was duly noted in some log book. From a bird’s eve view, a 
description could include movements of the forces geographically, and rates of change 
of both forces. In a global view, the overall evolution of the battle in time, and the 
search for underlying patterns in that evolution might be sufficient to describe the 
battle. 





Figure 3.1 Coarse Grid Representing Geographic Position. 


V7 


At each stage, or level of command, there is a need for a differing view of the 
battle, since at each level, the information needed is different. At the lowest level, 
concern might be for resupplies, 1.e., ammo, fuel, requests for transportation or support 
external to the lower organization. At a higher level, concern is for allocating the 
resources among competing requests and determining the priorities associated with 
those requests. At still a higher level, only the developing outcome of the battle might 
be relevant. 

We will attempt to describe an intermediate (“mesoscopic”) level between the 
upper (commander's or “macroscopic’) level and the lower (units or “microscopic”) 
level which incorporates all the relevant information a commander would need in order 
to make his decision. In defining this level, a set of order parameters needs to be 
developed. Order parameters represent this mesoscopic level and are a specific 
aggregation of the microscopic or state variables. For instance, in the tank example in 
the previous section, the state variables represent the detailed information about the 
tank, 1.e. velocity, position on the battlefield, training level of the crew, ammo supply, 
etc. A simple example of an order parameter in this case might just represent the 
number of tanks in a particular cell of the battlefield. An order parameter model 
would then develop the necessary interactions among cells, resupplv considerations, 
capability degradation, etc. We will see an application of this through the examples 
discussed later. 

As a first representation, both forces could simply be described as a certain 
number of personnel. We are then interested in describing this battle, given a set of 
initial conditions, in terms of force losses per unit time. This is equivalent to a 
Lanchester approach with noise, alluded to earlier. This is only one of several wavs to 
derive a stochastic Lanchester equation. This is referred to as a Langevin equation in 


the physics literature. This could be mathematically described as shown: 


dX/dt 


i ay ee CG ane 


dY/dt = P(X,Y) + 2 (X,Y) , (3.1) 


where X= the number of blue forces available to engage the Y forces, Y= the number 


of red forces available to engage the X forces, f *Y= functions relating average 
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numbers of blue forces and red forces, g*¥= functions multiplying (square-root) 
variances of the background noise, ai = coefficients relating rates of blue force and red 
force losses, and n= background noise. For example, functional forms might be 


ree 


ays ate airY : (3.2) 


fy = a5 1X aie a4 Y , (3.3) 


With similar equations for g*”. 

To attempt to solve this equation, we could put it on a computer, introduce some 
random noise (via a Monte Carlo simulation) and the aggregated output of many runs 
would give us at any time, t, via a probability distribution, the level of blue and red 
forces and any other variable which is dependent on these, such as force ratio, 
surviving maneuver force ratio, or some equivalent descriptive variable in which we are 
interested. 

We have selected the form (equation 3.1) because the current state-of-the-art 
mathematical physics then permits us to develop extremely general nonlinear means 
and variances into representations suitable for analysis by methods of statistical 


mechanics. 


C. EXAMPLES OF ORDER PARAMETERS 

To better understand the concept of order parameters, let us take a physical 
system as an example, a gas confined to a box in thermal equilibrium. The gas can 
obviously be defined at a microscopic level by representing it as a collection of atoms 
with certain velocities, relative positions, collision rate with other atoms and the walls 
of the box, and some internal energy state. In analogy to the battle, the atoms would 
represent the individual personnel, their velocities and positions corresponding to their 
movements and geographical positions on the battlefield, and the collision rate could 
correspond to the engagement rate with the enemy. The internal energy state could 
relate to the amount of ammo, firepower, and possibly training level of the individual 
combatants. However, on a more global level, there is a pressure associated with the 
gas, a temperature, and a volume. One of the problems with the order parameter 
concept is to find these global variables associated with combat and relate them to 


something of use to the commander. 
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The order parameter concept can be used to describe systems far from 
equilibrium. Combat is obviously a system far from equilibrium, except possibly in 
some isolated cases. Since the objective of both commanders is to accomplish some 
mutually exclusive mission, the system will tend towards a solution which favors one or 
the other commander. In analogy with our gas in the box, suppose we lower the 
temperature of the gas. At a certain temperature, the gas becomes a liquid which 1s 
called a phase transition to a long range collective order. The question is then: Is there 
an analogous “phase transition” associated with our forces and how do such collective 
patterns of information represent themselves? At what “temperatures” related to the 
size of the g*** functions does this occur? Is it a unique phenomenon, i.e. does it only 
occur at this “temperature”? What if we change the volume of the box, is there then 
some “phase transition” associated with our forces possibly relating to the change in 
geographic area of the battle? Are the order parameters of our physical system 
transformable to some similar order parameters of battle? 

In answering these questions, we can arrive at some understanding of combat 
and relate this to our understanding of other physical systems which undergo the same 
or similar transformations when the state of the system is changed. 

As a start, and following Bretnor [Ref. 6], two order parameters that seem likely 
are the destructive force and the vulnerability of the force. Destructive force is defined 
as the amount of combat potential which can be delivered to the enemy in order to 
destroy him. It includes the training, the readiness, the sustainability, the morale, the 
Weapons mix, etc. of the force. It 1s obvious that these factors do change during the 
course of battle, and that their level certainly would indicate the success or failure of 
combat. The vulnerability of a force are those factors which degrade the capability of 
the force, 1.e. position on the battlefield (terrain factors such as line of sight, cover, 
concealment, protective armor, etc.), lack of morale, discipline, or training, etc. As you 
can see, the vulnerability of a force 1s in some ways a degradation of the destructive 
force, yet they are not totally the opposite of the other. For example, a force in the 
open would have more vulnerability than one under cover, yet they would have the 
same destructive force. There are other examples, the point being that they are 
distinguishable order parameters, although we could effectively model the force using 
the destructive force and modifying it so that it is somewhat degraded when its 
vulnerability is increased. Nonlinear functions of the order parameters will be used to 
model scenarios in which the objective of the commander would be to attack the others 


vulnerability while exploiting his own destructive force. 
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IV. MATHEMATICAL FORMALISM OF THE MODEL 


A. INTRODUCTION 

In this chapter we develop the mathematical formalism of the model, and 
introduce mathematical objects such as the Lagrangian and the path integral. We will 
show there exist equivalent representations among the Langevin, Fokker-Planck and 
Path Integral descriptions of a stochastic system [Ref. 1,7]. We begin with a simple 
one order-parameter, non-linear model. The linear model is a special case of the non- 
linear model. We then fully develop the two order-parameter model which is used in 
the remainder of this thesis to illustrate the path integral method. Generalization to 
many order-parameters is made and included for completeness of the description. 
Assumptions of the model and their significance are given. The primary assumption 1s 
the requirement for a Gaussian-Markovian system. Finally, the relation to classical 
deterministic systems and the usefulness to classical statistical systems of the 


Lagrangian will be discussed. 


B. ONE ORDER-PARAMETER, NON-LINEAR MODEL 
The one order-parameter (1OP) model is not particularly useful in describing 
combat, but we present it for completeness and ease of illustration of the general 
model. The generalization to two or more order parameters can be easily made. Order 
parameters we could use are the ratio of the forces or the difference of the force levels. 
We begin by labeling our order parameter X, for example, the ratio of Blue forces 
to Red forces. We are interested in how X changes with time. Within a time 


increment, At, we can write 
A(t+ At) - X(t) = At f(X(t)) , (4.1) 


where f(X(t)) is some function to be fit. If At is assumed small and X is assumed to be 


continuous, we can then write 


X = dX/dt = f(X(t)) . (4.2) 
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In the context of describing combat equation 4.2 is referred to as a Lanchester 
equation and simply represents the mean or expected path of the order parameter for a 
large system. 

We now want to add a noise term to equation 4.2 so we can model the severely 


stochastic nature of combat. Hence, the change in X can be written 


X = X(t) + g(X(t)n (4.3) 


Where yj is the background noise with zero mean and variance=1 (assumed) and 
g°(X(t)) is the variance, which is not necessarily a constant. We also assume that n is 
Gaussian-Markovian (normally distributed “white noise”). The assumptions will be 
discussed in Section E. Equation 4.3 is referred to as a Langevin rate equation in the 
scientific literature, but we will refer to it as a generalized stochastic Lanchester (GSL) 
equation in the context of describing combat. This is only one way of obtaining a 
stochastic Lanchester equation. I.e., there is a many-to-one mapping of deterministic 
systems into stochastic contexts. 

Associated with this GSL is a Fokker-Planck equation [Ref. 8] defining a 
differential equation of the conditional probability distribution, P(X(t+ At)|X(t)), given 
as 


6P/dt = - G(fP)/OX + 1/267(g2P)/6X2 + VP. (4.4) 


The function f represents a drift term and g* is the diffusion term of the probability 
distribution P(X(t+ At)|X(t)).. Sometimes a potential term, V, is present, which 1s often 
useful to analytically enforce boundary conditions. | 

Another representation exists for describing P(X(t+ At)|X(t)) [Ref. 8]. For small 


time increments At, we assume 


P(X(t+ At)|X(t))= 1/(2mg2At)!/2 exp(-LAt) (4.5) 


Where L = (X - f)*/2¢7 is the Lagrangian of the system, i.e. a Gaussian form for this 
conditional probability, with mean (X(t)+fAt) and variance g7At, as follows for short 


times from equation 4.4 . It must be emphasized equation 4.5 is the short time 
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conditional probability distribution of P(X(t+At)|X(t)). With this representation, it is 
possible to obtain a long time conditional probability distribution P(X(t)[X(ty)) 
through a path integral description. This is defined as (more precisely defined in 
Langouche, et. al. [Ref. 7] and Schulman [Ref. 9]) 


P(X(t)[X(ty)) = foo dX ardXpoag’ | a At (4.6) 
x P(X(t)|X(t-At))P(X(t-At)|X(t-2At)) 


X++ +X P(X(ty + At)X(ty)) , 


= J-- JDXexp(- YAtL,) 


where 
DX = (2mg,7An 2 T(2mg_2Ary!/2 ax, 
n=1 


t,= tyt+nAt and t=t,)+sAt where t, are the intermediate time increments in the limits 
s— 0 and At —0. Equation 4.6 is called a path integral and is recognized as simply a 
Chapman-Kolmogorov equation. With the path integral, given some initial state at to, 
X(tg), We can determine the probability distribution of X at some later time t. The 
path integral is discussed in Appendix B. 

The purpose of the previous discussion was to show the relationship between the 
GSL equations and the path integral description of combat. This will also be shown 
for the two order-parameter and the many order parameter models, but in the actual 
development of the model all that is required is a functional form of the Lagrangian. 


This will be discussed further in section F. 


C. TWO ORDER PARAMETERS, NON-LINEAR MODEL 
We now develop the path integral representation of combat for two order 
parameters. We begin, as before, with a set of Langevin equations, show the related 


Fokker-Planck equation, and finally the path integral description. Our emphasis is on 
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the formulation and the notation of the path integral, whereas the Langevin equations 
are used to support our intuition. 

Suppose we are interested in our own force level and that of the enemy. We will 
use these as our order parameters and denote their level by X(t) and Y(t) for Blue and 
Red forces, respectively. 

As before in the 1OP model, we will be interested in the change of X, and Y with 


time according to 


X(t+ At) - X(t) = Atf X(t), Y(t) , 


Y(t+ At) - Y(t) = Atf2{X(t),(¥(t)] , (4.7) 


where the f! , i= 1,2 are some functions to be fit, and At is some small time increment. 
If we assume continuity of the order parameters and for small enough At, we can write 
equations 4.7 as 


X pl 


dX/dt 


Y 


dY/dt =f? , (4.8) 


These are the Lanchester equations (deterministic). 
e ® 
We now assume that multiplicative noise (Gaussian-Markovian on X and Y) is 


present and the order parameters are now modified according to 


Se eee 
Y=f2+ g? i+ 270° (4.9) 


where the ght are functions multiplying the variance of the background noise. If the 
it 
g 


the q .’s are assumed to be zero. We will also assume the number of noise terms 1s 


. were constants, then the n,’s would simply contribute “white noise”. The mean of 


equal to or greater than the number of order parameters [Ref. 10]. Equation 4.9 is our 


generalized stochastic Lanchester equation (GSL) for two order parameters. 


24 


The Fokker-Planck equation describing the evolution of the conditional 
probability distribution P(X(t+ At)|X(t)) where X(t+ At) = {X(t+ Ab), Y(t+At)} is 


OP/ét = VP + G&(-g# P)aM+ 1/26¢(gHYPy@MYSME , (4.10) 


-_ 2 — . . : 
where M! =X, M“=Y and V is a potential used to add constraints on the order 
parameters or to simulate boundary conditions. The indices ft,v = 1,..., N where N 
is the number of order parameters (2 in this model). The git and gHY are different 


functions from the git in the GSL and are defined as 


gt = f+ i/2g¥ dgh GM” , (4.11) 


giY = gh ov . (4.12) 


We are now using the Einstein summation convention, whereby repeated indices in any 
term imply summation over those indices. In the 2O0P model these are, for example 


when w= 1, 


. Gg! dg! 6g! Gg! 
Los ff + yagt); =++1/2¢! 24 12g?) —++1/229°, — . 4.13 
= Ge ee gee or ay ie 





(The compactness of the Einstein summation convention is evident here) and for p= 1, 


v= 2, 
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meee) 2 lez 
g = 2,8) +238. - (4.14) 
Note the ght which are the variances of the microscopic noise sources are summed over 
and contribution from individual sources need not be fit in the path integral 
description. The path integral description of equation 4.9 is 
$ 


P(X(t)|X(tp) = f° * ‘[DXexp(- YAL,) 
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/ sfio™ = f> 2N -1/ 
DM = g,!/* (20 At)? qT gt? Ona M2 ME 


fame + YAMF, ; MB, = ME, MEL) = ME, 


Pet 200 eee. ie ane 
Zayem (') (4.15) 


eae cl (2nv)n 


This is the long time conditional probability distribution of our order parameters. The 


short time conditional distribution is 
P(X(t+ At)[X(t)) = g!/2 (2mAty!/2 exp(-LAt) (4.16) 


where L and g are as defined above. 


This description is correct as long as we adopt an Ito or pre-point discretization 
of our order parameter, 1.e. 


MFP (t.) = ME, 
(4.17) 
MM (t,) = (MP MEM ty) 
and t, = tg tnAt. This affords us the luxury of a relatively simple Lagrangian. 


There exists a mid-point or Stratonovich discretization of the order parameter given by 


MA(t_) = 1/2MF + ME ) : 
(4.18) 
MF (oa (MP eee Mia ae 


This induces a curved or Riemannian space on the order parameters with the 
subsequent requirement of additional terms being added to the Lagrangian. The 
presence of the noise actually induces the non-Euclidean geometry of the fl-space and 


the variance gH” is the inverse of the ft-space metric, Envy: The benefit of having a 
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mid-point discretized Lagrangian is that the associated Euler-Lagrange equations 
determine a variational principle. This allows us to derive a most likely path of the 
order parameter, without doing a full calculation of the long-time probability 
distribution [Ref. 1,7,11]. 

The preceding discussion was not meant to be rigorous, but to point out the 
subtleties in actually evaluating any of the functional forms. For simplicity we will use 
the pre-point discretization and not carry any of the Riemanian terms. An example of 
the two order-parameter model with an explicit form for the Lagrangian is given in 
Chapter 6. 

Although we have developed a Lagrangian from the GSL equations in the 
preceding discussion, it is not necessary to do this in general, i.e. we can begin our 
model by assuming a functional form for the Lagrangian without having first written 
down the associated GSL equations. This is an important point. There is an algebraic 
relationship between the Lagrangian representation and the GSL representation (under 
the assumptions) and one could, in principle, derive the GSL from the Lagrangian and 
vice versa. There exists a large bodv of literature on combat modeling with Lanchester 
equations and thus the experience gained using that approach can be transformed to 
the Lagrangian approach. There also exists a large body of literature dealing with the 
applications of the Lagrangian approach to other large, complex, physical systems 
which can then be directly used to provide physical insight into the problems of 


modeling combat. 


D. MANY ORDER PARAMETERS, NON-LINEAR MODEL 

The extension to many variables can be made [Ref. 1,12]. Suppose now we are 
interested in modeling the spatial-temporal patterns of the order parameters and not 
simply the temporal patterns as before. To extend the 2OP, where ht = 1,2 , we now 
let w = 1,...,N, where N is the number of order parameters we want to model. For 
one example of a many parameter model, suppose we divide the battlefield into distinct 
cells, labeled by @= 1,...,m. For example, in each cell we would examine the Blue 
and Red force levels as composed of tanks and personnel and as shown in Figure 4.1. 
The fH@ where ya forms an enlarged index of the 'X@ variable-space, and V can 
now incorporate NN (nearest-neighbor) interactions and N2N (next-nearest-neighbor) 
interactions to account for external forces, such as higher level constraints, resupplies 
from adjacent units, actual movement of forces from cell to cell, etc. The model would 


be as follows, 


Za 


P(X(t)IX(tg) = fr “fOXexp(- PAL) , 


i N m 
DM = g,!/2 (aman? fe V2 EFT (2mAty 2? aM 
uw ¢ 


n=1 


g, = det (fuy op) (4.19) 
L = 1/2(Mbe pH ey poy? _ £ VB) (4.20) 
M = {M/® jp=1,...N a=1,..4m n=1,.. .,5} (4.21) 


It must be emphasized that we are only assuming a Gaussian distribution of the rate of 
change of the variables in time, and that the spatial distribution could be non- 
Gaussian. This is a modeling consideration when deciding on a functional form of the 
Lagrangian. It should also be emphasized the distribution is only Gaussian in the 
short time, and only in the post-point value of the variables, whereas the long time 


distribution could be any distribution. 


Ek. ASSUMPTIONS 

The primary assumption of the general model is that the system to be modeled is 
a Gaussian- Markovian system in the rate of change of the variables. This means there 
must be sufficient order parameters available to describe the system as Markovian, 1.e. 
that the future state of the system only depends on the present state. This assumption 
comes into play in describing the short time conditional probability distribution of the 
order parameter. The Gaussian assumption states the short time conditional 
probability distribution is Gaussian in the post-point variables. This is a standard 
assumption made when dealing with many stochastic models. It is hoped the 
Lagrangian or path integral representation is very robust with respect to these 
assumptions. This means if the noise is non-Gaussian or there are not enough order 
parameters to ensure a Markovian description, then we can still obtain a reasonable 


path integral representation with these assumptions. 
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Figure 4.1 Battle Grid Showing Blue and Red Force Parameters. 


FL. INTERPRETATION OF THE LAGRANGIAN 

Suppose We have a functional form of the Lagrangian which can be plotted as in 
Figure 4.2. We now show how the Lagrangian contains information about the svstem: 
most likely states of the system: a measure of the risk associated with that state; and a 
measure of the transition probabilitv between most likely states. 

The minima of the Lagrangian correspond to the most likely states of the svstem. 
We assume we are looking only at the short time probability distribution, eae 
Lagrangian contains a widelv varying expression containing factors of the difference 
between the actual state and the average or mean state, te. L © [Np4 yy - (A, + fAt)}- 
where the term in the parenthesis is the past state corrected for the drift. Therefore, if 
the actual state is much different from the average state. then L will be a relatively 
large value compared to the other terms, and the corresponding value of the 


probability distribution will be correspondingly exponentially small. 
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Figure 4.2 Plot of Lagrangian vs Order parameter. 


The minima are shown as points A,B and C in Figure 4.2 and represent the most 
likely states of the system) at time t. 

The Lagrangian also contains the uy term Which is the metric of the order 
parameter space, 1.e. it is a measure of the distance in this space. This space is curved 
for all metrics which are not constants, and we then have a short time Gaussian 
distribution in curved space. 

The ny term is also related to the variance gltY = (Buy) , of the underlying 
mucroscopic sources of noise which is a measure of the “width” of the minima. The 
minima width of a most likely state can be seen as a “degree of risk” measure 
associated with that state, i.e. the wider the minima the larger the confidence interval 1s 
for that state. For example, look at the minima at point A. If we were trying to 


forecast the value of the order parameter X, then we could only say it is betwcen X= 9 
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and X=13, with any degree of confidence (to be defined). At point B we have a 
minima which is more sharply defined. We have most likely states of X=15, with a 
confidence interval of X=(14,16) and the g#Y would be a much smaller term. The 
smaller the g4” the sharper the width of the minima. A measure of the transition 
probability between two most likely states is the relative height of the “peak” between 
the two “valleys” of the minima. For example, point D is much smaller than point E, 
and if the state of the system was at point B, then a transition to state A would be 
more likely than a transition to state C. 

Granted the above discussion is very qualitative, but quantitative results can be 
obtained. However, a graphical device portraying the information contained in Figure 


4.2 would be more useful as a qualitative tool than a quantitative one. 
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V. DESCRIPTION OF THE METHOD 


INTRODUCTION 


We will now present, in outline form, the generalized procedure for obtaining a 


path integral representation of combat systems. 


Select, derive or develop the order parameters of the system thereby defining the 
independent variables 


Obtain sufficient empirical data from the system you wish to model 

Functional forms of the independent variables are developed in terms of 
theoretical parameters/coeflicients to be fit to data, to model means and 
variances 


Perform a_maximum likelihood fit of the short time probability distribution, 
fitting coefficients of the functional form 


Using the path integral technique, a probability distribution of the order 
parameters is found for long times. 


Perform sensitivity analysis 

With the probability distribution, you can then use the method to 

- Analvze budget decisions in terms of hardware/software purchases 
- Perform combat analysis for use in battle management 

- Determine the effect of proposed doctrinal changes 

- Perform “What if” scenarios for use in combat planning 


- Have additional input to your decision making cycle 


The method is an iterative process. We will collect some data, look for order 


parameters, attempt a fit, and if not very successful, try a new functional form of the 


Lagrangian. This process will continue until we have decided our assumptions are 


satisfied and we have a reasonably good fit to the data. Of course, after examining the 


structure of the Lagrangian and discovering that the model gives results which are not 


correct, then we must go back to the beginning or try a different functional form for 


the Lagrangian. We will now cover each of the steps in detail. 


B. 


ASSUMPTIONS OF THE METHOD 


Before we begin describing the method in detail, 1t would be appropriate to look 


at the assumptions. These assumptions will guide our development, and selection of 


the order parameters. Selection of the order parameters will then guide our data 


collection efforts. 


As stated before in Chapter 4, our primary assumption is that the system to be 
modeled is Gaussian-Markovian. This assumption was necessary in developing the 
path integral. This assumption has further consequences in defining the amount of 
data required in order to sufficiently model our complex combat system. 

There are several items which need to be addressed: 

e number of elements in the battle; these can be personnel, vehicles, aircraft, 
ships, etc. If they are to be used as an order parameter, then there must be 
enough individual elements to ensure approximate continuity. 

e number of runs of the experiment/war game/simulation; this is required in order 
to provide sufficient statistics for a good fit of the Lagrangian to the data, 
mainly in estimating the parameters of the g#, i.e. the variance, and the means 
gh. 

e number of order parameters; this is required to ensure the system is Markovian. 
If not enough order parameters are used, then even if the true system is 
Markovian, our model of it may not be Markovian, which could result in a bad 
fit of the Lagrangian. If the model is robust, then a good fit could still be 
obtained if the number of order parameters used is not too different from the 
actual number of underlving order parameters. There is a subtle but important 
feature of our modeling which helps to create a robust fit: when care is taken 
to handle all nonlinearities, e.g., including Riemannian terms in the midpoint 
discretized Lagrangian, equation 4.18, then the probability distribution is 
invariant under nonlinear transformations of the variables. (This is what 
induces the Riemannian geometrv [Ref. 10].) Thus we are really fitting a wide 
class of functional forms whenever we do one generic fit to the data. 

° “uncertainty principle’ O<t= At < I/L, L ~ <(Ax)>?/ (<(Ax)> At) 
where ~ means on the order of. This places a requirement on the amount of 
change in the order parameter in a particular time increment. This is also 
necessary to calculate the path integral. This states that in a time increment f, 
the average drift of the order parameter must be less than (or on the order of) 
the variance [Ref. 13]. Obviously, when considering actual data, Le. from 
operational exercises or combat, then At will become part of the data and then 


we can only do the best fit we can with the data available. 
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e Aggregation (co-location) of capabilities; it is possible this tvpe of model may 
not be adequate if the system is composed of only a few large distinct entities 
which have several capabilities. For example, a naval battle group may have 
large numbers of personnel and aircraft, but they are aggregated into a 
relatively few number of ships. Destruction of one ship’s capability may have a 
large impact on the outcome of the battle. It still may be possible, however, to 
model some aspect of the naval battle, for example the group’s outer air battle, 
which has sufficient numbers of elements 1.e. aircraft. There is more interesting 
work that could be done here but is beyond the scope of this thesis. 

In summary, we need a combat system which has a large number of elements to 
ensure approximate continuity, a sufficient number of order parameters, and a 
sufficient number of experiments to provide for a good fit of the Lagrangian. Once a 
Lagrangian is developed and coupled with the path integral, we will have a 
“propagator” to describe the time evolution of the system from any initial time t, to 


anv final time, t. 


C. SELECTION OF ORDER PARAMETERS 

The development of the order parameters is first dependent upon the system you 
wish to study. For example, if you were interested in the length of a battle as defined 
by some cutoff strength for either side, then the order parameters used nught be just 
the force levels of each side. If you were interested in the relation of C? to the 
outcome of a battle, then you might use some MOE of C? of either side, together with 
the force levels. 

Second, the order parameters used must satisfy the aforementioned assumptions. 
Obviously, an order parameter must be something which changes value during combat 
and cannot be a constant or a very slowly changing variable. Otherwise there would 
be no need to model that particular order parameter. 

Some examples of order parameters (per side) are: 
number of vehicles (tanks, trucks, aircraft, etc.) 
number of aircraft 
number of elements firing (artillery/tanks/aircraft) 
number of shots fired 
number of tanks (armor study) 


supply/logistics availability 


ee ee 


troop carrying capacity (helicopters/tactical troop carriers) 
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8. geographical position on the battlefield (appropriate for motorized infantry) 
9. number of aircraft/sector ( outer naval air battle) 
10. log(bits) of information used in communications 
Through the appropriate use of the order parameters and the path integral other 


questions such as the outcome of the battle and duration of combat can be answered. 


D. DATA COLLECTION 

The first step in the analysis is to obtain empirical data. This could be the result 
of simulations, war game results, field exercise data or data from actual combat. Each 
of the data sources has its advantages and disadvantages. If the data were available, 
the best source would be from actual combat. Obviously, there have been no large 
scale wars recently, but that does not prevent us from doing analysis on past wars. 
This will not be the approach here. There is data available from field exercises, but it 
is not in a suitable form for analysis at the present time. This includes data from 
CAX’‘s (Combined Arms Exercises) which are held at Twentynine Palms, California by 
the Marines, or exercises conducted at Fort Irwin by the Army on their calibrated 
range. This could be done at a later ttme. War games are the next best place to 
obtain data. Several War Games available at NPS are JANUS (after the mythological 
two-faced god) and IBGTT (Interactive Battle Group Tactical Trainer). TWSEAS 
(Tactical Warfare Simulation, Evaluation and Analysis System) is a war game available 
at the Marine Corps Development Center in Quantico, Virginia. These simulations are 
discussed in the Appendix C. Other simulations available are CARMONETTE, 
SOTACA, and FOURCE to name a few. A brief description of each is also included in 
Appendix C. 

For the purposes of constructing a statistical mechanics model of combat, many 
trajectories of the order parameters are needed. What do we mean by trajectories of 
the order parameters? In the space defined by the order parameters, a point represents 
one possible state of the order parameters. A path which connects the initial state of 
the system to some final state is called a trajectory. The trajectory represents one 
possible realization of combat. For a good fit of the model to the data many 
trajectories, and therefore, many stochastic experiments are needed. This will naturally 
leads us to select a war game or simulation as a source of data. In these cases, many 


experiments can be completed with variations in the noise. 
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E. DEVELOPMENT OF THE LAGRANGIAN 

The next step in the analvsis is to determine a functional form of the Lagrangian. 
The simplest and most versatile form is a ratio of polynomials. This defines a Pade 
approximant form and is suitable for approximating many functional forms. The 


Lagrangian 1s 
L = 1/27MH. fe apy” -£%) (5.1) 


where we need to assume functional forms for the fF , w= 1,...,N (Nis the number 
of order parameters) and the variance g4Y . Since Zyy is a metric it must be positive 
delinite, -emdet (Syy)> QO. Except for this one condition, the Lagrangian can be of any 
form. Obviously, we would like to keep the form simple, yet model the data 
accurately. This might require several iterations of: 

1. select a functional form of the Lagrangian 

2. performing a maximum likelihood fit to determine the unknown coefficients, 

3. testing the fit of the Lagrangian with the data 

4. andif not satisfactory, go back to 1. | 
If, after several iterations, a good fit has not been attained, then we must look at our 
data to ensure we have satisfied our assumptions. This could be one test to see if we 


satisfy our assumptions. 


F. MAXIMUM LIKELIHOOD FIT 

After deciding on a set of order parameters defining independent vanables and a 
functional form of the Lagrangian, we are now in a position to estimate the 
parameters/coefficients of our Lagrangian. This will be accomplished by using a 
maximum likelihood fit. 


The short time conditional likelihood function, M, is defined to be 


(2nAt)}/2 gi/2 exp(-LAt) = P(X(t+ At)[X(t)) , (5.2) 


= 
ll 


te 
| 


= 1/~ME. FB epy(M _fV) (5.3) 
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where the other variables are as previously defined. Our data is in the form of I runs 
for JAt time periods, where one run ‘is one realization of combat and J is the duration 


of combat. Our maximum likelihood function for many runs now becomes 


Ea 
i 
0 et 


a | j 
IT (2m At)! /? gl’? exp(-L..At) 
l j=l j IJ 


T 
= LIP(X;(ty + JA0IX, (ty + G-DAt) x oa P(X,(ty + At)|X.(ty)) (5.4) 
where X(tp+jAt) = {X(t)+jAt),Y(ty+jAt)} , 


and where the data is used to calculate specific values of the X’s. We now wish to find 
the parameters in the Lagrangian which maximizes this function. To do this we first 
take the logarithm of M’ to accommodate computer requirements on acceptable ranges 
of numbers which can be processed. Since the logarithm is a positive monotone 
function, the maximum of In M’ will be in the same location as that of M’. Therefore we 


now wish to maximize 


T J 
In M = me p> [-1/2ln(27At) + 1/2In g.- L,At J (66) 
T J 


=. a ‘. [-1/2In(2) +1/2In(At) - 1/2In g, + L, At] 


The maximization of In M’ is equivalent to a minimization of -In M’. Also the constant 
1/2ln(27) can be deleted from In M’ since it will not affect the location of the nunimum. 
In general, At will be part of the data and thus we will not drop this term. Therefore 


our problem is to locate the minimum (global, if possible) of 
N= )Y [1/2ln(At) + L,At - 1/2ln g:, ] (5.6) 
Current algorithms for solving non-linear minimization problems are 


deterministic and only guarantee local minima. Therefore we have developed a version 


of the simulated annealing algorithm which guarantees convergence to the global 
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minima if certain conditions are satisfied. The algorithm uses a random sampling rule 
to select points from the parameter/coefficient space, and a criteria on which to base 
acceptance of that point as the new state of the fit. We use a Cauchy distribution tied 
into a “temperature” on which to sample the space, and as the temperature is lowered, 
the search is more localized. The Cauchy distribution (a long tailed, © variance 
distribution) is used so that the localized search does not get trapped in a local minima, 
and there is some probability of leaving to find better optima. The Boltzmann 
distribution (from the Metropolis algorithm) is used as a criteria on which to accept 
the point. This distribution is particularly useful in our case, because the cost function 
is in the form of a Boltzmann distribution and thus there is a more efficient mapping to 
the parameter space. A description of the algorithm, the necessary conditions, and a 
FORTRAN program of the simulated annealing algorithm are contained in Appendix 
A. 

There are a few subtle points to mention concerning data, variables, parameters 
and spatial dimensions. The paths which are generated from a source of data are 
equivalent to a set of likely paths which can be sampled from the theoretical 
distribution of the variables. This is usually what is meant when discussing sampling 
statistics. The data are a sample from some theoretical distribution which is at present 
unknown. The parameters are the coefficients in the functional form of the Lagrangian 
and the Envy metric and are to be estimated from the data. At the time we are fitting 
the data, M becomes a function of the parameters with the variables being the data. 
Once the Lagrangian is fit to some set of the data, it then becomes the fitted 
distribution of the underlying order-parameter variables which attempts to mimic the 
unknown underlying theoretical distribution. The “dimension” of the space of variables 
is determined by the number of order parameters being considered. If two or three 
spatial dimensions are being considered, then typically separate cells within this space 
would contain independent order parameters, which would be functionally tied to the 
order parameters in other cells by the functional forms used to define the multivariate 
drifts and diffusions. 

It should be stressed that fitting the Lagrangian does not mean separately fitting 
the means and diffusions to the same accuracy. That is, we are fitting the functional 
form of the Lagrangian to the data and NOT finding a set of parameters/coefficients in 
the drift and diffusion terms of the data. Thus there may be a “wide” discrepancy 
between parameters/coefficients using similar data, but it is the resultant fit of the 


Lagrangian to the data that is most important. 
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G. PATH INTEGRAL TECHNIQUE 
The most important part of the analysis is now complete, i. e. developing a 
functional form of the Lagrangian. With the Lagrangian in hand we are in a position 


to calculate the long time conditional probability distribution P! defined as 
P! = P(X(t)|X(tp) = f-°-f DX exp(-LAt) , (5.7) 


where ty) < t and t can be any time chosen in the future. This is the path integral and 
acts as a propagator of the system, i.e. once the path integral is known, given any 
initial state, the probability distribution of the order parameters at any other time can 
be calculated. For simple systems, Monte Carlo techniques are used for multi- 
dimensional systems. However, there is only one method, recently developed, that has 
proved to be accurate for a wide range of nonlinear, nonstationary problems, such as 
those we expect to be present in combat systems. At the present time, using this 
method, the path integral has been calculated in one dimension [Ref. 13,14,15] for 
many highly nonlinear systems, and it has recently been extended to two dimensions 
[Ref. 16]. Work is ongoing at Lawrence Livermore and Sandia National Labs to 
develop an algorithm for the many dimensional case. However, meaningful results can 
still be obtained by examining the static Lagrangian where X ~ 0. This gives some 


indication of the characteristics of the system before doing long calculations. 


H. VALIDATION AND SENSITIVITY ANALYSIS 

Now we are in a position to check the sensitivity of the model to changes in the 
data. The Lagrangian essentially contains all the elements of the model. This 1s why it 
is SO important to obtain a good fit. This analysis could be done in two ways: 

1. With many runs of the data, we can separate (randomly) a set of runs to do our 
fit and then validate the fit using the remaining set of runs. It is not clear at 
the present time how to determine a good fit but a method that seems 
reasonable is determining deviations from the most likely path-in the following 
manner. First, at each time increment, calculate the sample variance of the 
data. This will be our weighting factor. Next, calculate the distance between 
each of the data points and the most likely state calculated from the 


Lagrangian. This will be our “width”. We then form the following statistic: 
1(n,(n,-1)) L(Aw,)7/s7, = (nae LOR)", (5.8) 
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Where 4 = Li(n,- DY (%,-x,)* x, is the most likely value at time t, x, is the 
sample mean at time t, n 1s the total number of data points, n : is the number of 
experiments and n, is the number of data points time t. The distributional 
properties under our assumptions of this statistic remain to be seen. One 
property which is evident is if the long time conditional distribution is 
symmetric then x, and x, are equal and our statistic will be unity. Another is if 
the data were to fall within one standard deviation s of the most likely path for 
all time, then the statistic will be close to 1. In this case a good fit would be 


evident if the statistic were “close” to l. 


to 


Using a set of runs of a particular scenario, fit a Lagrangian. Now change a 
parameter of the scenario, for example, starting force levels. Next use a 
goodness of fit test such as described above to check the sensitivitv of the 
Lagrangian to the change in the scenario parameter. This could be done for 
several parameters or several iterations of the same parameter. 

Only after a full sensitivity analysis has been completed and the Lagrangian can 


be shown to be robust, can we say the Lagrangian accurately models the scenario(s). 


I. OPTIONAL USES 

We now present several generic uses of the path integral method and show some 
of its usefulness and versatility as a combat model. 

1. Procurement Decisions 

Suppose you are a commander in charge of procurement decisions. You are 
faced daily with comparing Cc systems to tanks to aircraft to field artillery pieces and 
up to now have relied mainly on subjective or qualitative models. With PIACA, the 
comparison between different items of equipment/weapons is summarized into 
comparing probability distributions and actually compare the items value in changing 
the outcome of combat. 

For example, suppose we are interested in comparing the relative worth of 
purchasing a new C system or purchasing more tanks. A study could be conducted as 
follows: 

|. Obtain a sufficient physical model of the C? system which can be simulated or 
used in a war game. A similar requirement exists for the tanks. Although this 
seems somewhat complicated, this step is usually completed for most 


comparisons. 
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2. Once a sufficient physical model has been obtained, conduct many runs of a 
simulation/war game using first one system, then the other using various 
appropriate scenarios. 

3. We now fit a Lagrangian to both data sets and develop a short time conditional 
distribution and a long time propagator or path integral. 

4. We now have common objects, the Lagrangian and the path integral in which 
to compare the two different systems. We might then compare the most likely 
states of the two Lagrangians and determine if they are satisfactory. Or we 
might compare the long time distributions or develop some MOE which 
combines features of the Lagrangian and the path integral. 

The are several advantages to this approach. First, fewer simulations are 
required since the information present in the scenario is captured by the Lagrangian 
within a relatively small number of runs. Second, once the fit has been obtained it is 
possible to perform additional sensitivity studies using the Lagrangian without 
requiring more, expensive simulations. This leads to a quantitative and objective tool 
which can be used by procurement managers. 

2. “What If” Scenarios in Combat Planning | 

Now that we have seen how to use the method in procurement decisions, it 1s 
not much of a jump to the use in combat planning. There is a common thread to the 
method, the comparison of similar quantities, the Lagrangian and the path integral. 
For illustrative purposes we present an example of the method in combat planning. 

We now suppose we want to examine the effectiveness of several combat 
plans. These are not as dissimilar objects as before, but we now have an additional 
objective evaluation which we may use. A study would follow similar lines as before: 

1. Perform many runs of simulations for each combat plan. This alone would aid 
in understanding the significance or utility of each plan. 

2. Fit a Lagrangian to each data set. Perform any required/desired sensitivity 


analysis. 


od 


Calculate the long time distribution via the path integral. 
4. Develop any MOE’s or compare the distributions in a qualitative manner. 
Another advantage of the method is the user becomes more involved and 1s 
able to see the consequences of each plan more fully. Examining the structure of the 
Lagrangian enables him to see transition points in the developing outcome of the battle 
which might be critical points. With this information, he can then make more 


informed decisions regarding the evolution of the battle. 
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3. Doctrinal Evaluations 

This use again follows the same procedure and we will suggest several 
examples. The procedure followed is the same as above. 

One example could be to examine the effect of a change in armor tactics due 
to a small change in the weapon system. A large change in the weapon system might 
require a new study and therefore a new Lagrangian. 

Another example is to look at a change in defensive tactics such as the 
difference between using a line defense and using a dispersed defense. 

We could look at a change in helicopter tactics as a final example. 

I have attempted to give the flavor of the possibilities of the method. It is 


extremely rich in applications and much interesting work can be accomplished. 


J. SUMMARY 
We have provided a general methodology for developing a statistical mechanics 
model of combat in this chapter. It is as follows: 
e Derive or develop the order parameters of the system 
e Obtain sufficient empirical data from the system you wish to model 


e Functional forms of the order parameters are developed to model means and 
variances 


e Perform a_maximum likelihood fit of the short time probability distribution, 
fitting coefficients of the functional form 


e Using the path integral technique, a probability distribution of the order 
parameters is found for long times. 


We have suggested several possibilities of the method for use in procurement, 
combat planning, and doctrinal evaluations. The utility of the method is being able to 
extract the essence of a combat system and representing it by a functional form, 1.e. the 
Lagrangian. Once the Lagrangian has been found, it then becomes possible to predict 
future outcomes with some degree of statistical (un)certainty. 

We have also argued that although the Lagrangian is specific for a particular 
scenario, it may be robust enough to small changes in the scenario. A collection of 
Lagrangians is possible for example, to describe differing scenarios such as cold 
weather, desert, or jungle environments, defensive versus offensive tactics/strategies, 
differing force structures and differing advantage factors, 1.e. either for or against the 
enemy. This collection after suitable testing can then be incorporated into the decision 
aid PIACA, for use by the commander and his staff. 
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VI. THREE EXAMPLES OF THE METHOD 


A. INTRODUCTION 

In this chapter we look at three examples of the use of the path integral or 
Lagrangian representation of combat. The first is a one order-parameter model with 
constant variance. This will lead to a quadratic Lagrangian which will imply a 
Gaussian distribution at each time step with means and variances being functions of 
the order parameters. The second example is a two order-parameter model with 
multiplicative noise. The fit becomes more difficult but results can be obtained. For 
both of these examples, simulated data from the associated GSL is used and it will be 
shown that the coefficients used in the GSL can be obtained from the fit using the 
Lagrangian representation. These two examples are provided to illustrate the 
relationship between the GSL and the Lagrangian representation, and as a test of the 
computer code of the maximum likelihood fit program. It must be emphasized that we 
will be fitting a Lagrangian to the data and not merely to find coefficients which are 
exactly the same as those in the GSL. Thus we will compare the theoretical Lagrangian 
(derived from the GSL) and the empirical Lagrangian which is fit from the data. In 
the third example we show how to proceed with the method when given a set of 
empirical data which was taken from the war game JANUS. Due to time constraints it 


was not possible to perform the necessary fit to the data. 


B. ONE ORDER-PARAMETER MODEL 
Although as stated before, the one order-parameter model may not be a good 
model for combat, we present it here as a simple exposition of the mathematics and the 
methodology. 
1. Data Collection and Order Parameter Used 
The data was generated from a generalized stochastic Lanchester equation of 


the form 
SaaX + en. (6.1) 


The n are distributed N(0,1). (Note: This is, in fact, an Ornstein-Uhlenbeck process.) 
The APL program used to generate the data for this example is given in Appendix D. 
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The order parameter represents, for our example, the Blue force level, which 1s 
being attrited through some process such as an essentially constant Red force level, 1.e. 
a = a@yY. The data consists of 20 runs of the simulation for total time increments of 
s= 20, and At = 0.1. Therefore the total time of one simulation represents, sAt, or say 
2 minutes. We are assuming there is no explicit time dependence in the Lagrangian, 


and any time dependence is implicit in the variable X. Our GSL for the data is then 

X = -1X + 4 where n ~ N(O,1) . (6.2) 
For generating data we take the form 

X(t+(j+ 1At) = X(t + jAt)+(-O0.EX(t + jAt)+n)At (6.3) 
Sample trajectories for 20 runs are shown in Figure 6.1. It should be noted here that 
all figures in this this chapter used the powerful statistical and graphics capability of 


GRAFSTAT. Figure 6.2 plots the sample distribution for each time increment and 


also connects the means of the distributions. It is obvious here that the expected value 


path is the solution of the deterministic equation, te. <X-aX = n> = <X> - 
® ° 
ax<X> = 0 — <X> = a<X> since <n>= 0. The notation <argument> 


signifies the expected value of the argument. This will be seen to be the case when the 
Lagrangian is quadratic. 

We now pick the short time conditional distribution representation to fit this 
data: 


P(X(t+ At)|X(t)) = 1/(2mg*At)!/? exp(-LAt) (6.4) 


Where L = (X 2 ha2e- = (X-a,X)/24,7 and where the a.s are parameters to be 
estimated. For this example, 


P(X(t+ At)/X(t)) = (2WAt)!/2 expf-1/2[X(t+ At)-X(t) 


(6.5) 
-(A, X(t))At}?/a,7At} 


se 


SSS 





Figure 6.1 Trajectory of Order Parameter X. 


and we can see this is in a Gaussian form with mean ft = X+(a,X)At , and variance 
a7 At . We also see the Lagrangian in equation 6.5 is in quadratic form. 
2. Maximum Likelihood Fit of the Lagrangian 
We will use a maximum likelihood fit to estiniate the parameters a, and 4, in 


the Lagrangian. Therefore we want to maxinuze 


= 
lt 
oo = 


1 


HL P(X(t+(j+ AVIX(t + jAQ) (6.6) 
j=1 : 


or equivalently to maxinuze In M, 


nM = Yo PY in P(X(t+(j+ YAIX(t + jAt) 
a J 
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DISTRIBUTION OF ORDER PARAMETER 





Figure 6.2 Distribution of Order Parameter X. 


= YY (/ain 2ndt) + 1/2in g? + L,At) 
1 J 


To use the simulated annealing algorithm described in Appendix A, we need our cost 


function (likelihood) in terms of a minimization. Therefore, we want to minimize 


Q -InM’ = ¥) Y' (1/2 In (2mAt) + 1/2In g*+ L;,At) 


SY (1/ain (2wAt) + 1/2ln aj? + 1/2AYX(t+ (i+ NA) 


X(t + jAt) -(a,X(t + jAt))Aq]} 


There is an algebraic relationship between the GSL and the Lagrangian. The 


coefficients/ parameters in the Lagrangian correspond to the coefficients in the GSL, if 
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defined properly. Therefore, a good fit would be obtained if the estimated parameters 
of the Lagrangian were “close” to the coefficients used to generate the GSL. However, 
it must be stressed again that we are fitting a Lagrangian to the data and not merely 
attempting to match the coefficients. We will compare the theoretical Lagrangian 
(derived from the GSL) to the Lagrangian fit from the data. 
3. Results of the Fit 

In addition to the above example, two other fits to data for different GSL’s 
were performed to truly test the simulated annealing computer code. The results are 
listed in Table 1. In defining the model we will use the terms linear or non-linear to 
describe the functional form of the drift and additive or multiplicative to describe the 
form of the diffusion or noise terms. As is evident in the table, the fitted coefficients 
were “close” to the generated coefficients. The true test however, was how well the 


generated Lagrangian and the fitted Lagrangian agreed. 


TABLE 1 
ONE ORDER PARAMETER RESULTS 


Model GSL Generated Fitted 
Coefficients Coefficients 
I Linear Additive X=aX+gn a =-0.008 A, = -0.0121649 
e=0.2 4, = 0.104862 
IT Non-linear X=aX+bX?+gn a=-0.008 4, = -0.0385419 
Additive b=-0.001 d, = -0.000818763 
g=0.2 4, = 0.104391 
HI Linear X=aX+gXy a= -0.008 4, = -0.0121295 
Multiplicative 2=0.01 4, = 0.00482471 


Before comparing, the Lagrangians had to be renormalized. Due to the 
underlying noise the Lagrangians could only be fit within an arbitrary constant. This 
constant was chosen so that the renormalized Lagrangians were unitless (this was 
arbitrary). A first choice was to use as a constant the Lagrangian evaluated at Ly 
(X=0, X= X_) where X, is the value of X which L is a global minimum. However, 


Ly(0,X%,) will more than likely be close to zero and this could cause numerical and 
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mathematical difficulties. Therefore, we chose to evaluate Ly, = L(0,X;) where Xy 
1. Re 2 > : : 

=e eee (X_)/ t. g°(X_) is the variance evaluated at X,- This seemed appropriate 
since we would like to incorporate into the normalization factor a measure of the 
curvature of the order parameter space and yet be “close” to the global minimum. 


Therefore 
Lp = L(X,X)/Ly(O,Xz) . (6.7) 


is our renormalized Lagrangian. These are plotted in Figures 6.3, 6.4, and 6.5. It can 
be seen that the agreement in the drift terms measured by the location of the minimum, 
are very good. Ihe agreement in the variance measured by the shallowness of the well 
was fairly good. This was for 20 runs of the simulation from time 0 to time 2 at 
increments of 0.1. The estimated coefficients were determined by having the simulated 
annealing program run for 10000 generated points and selecting the minimum value 
obtained. More work is being done here in running the program for 1,000,000 points 
and also performing the fit with more data, i.e. more runs and more time increments to 
see the effect of these modifications on the location of the minimum. 

These simple examples were meant to illustrate the principles of the method 
and thus simple Lagrangians were obtained, 1.e. those with only one minima. 
However, more complicated Lagrangians are admissible and this method is only limited 


by the ingenuity of the modeler. 


C. TWO ORDER-PARAMETER EXAMPLE 

The two order-parameter example is now given. The Lagrangian for this 
example involves a ratio of polynonmuals because of multiplicative noise. We again 
show the coefficients used to generate the data can be estimated quite well using the 
two order-parameter Lagrangian. This example tests our computer code, so that if the 
functions used in the Lagrangian are derived in a special way from the Langevin 
equations, then the short time conditional probability distribution will correspond 
directly to the probability distribution of the variables generated by the Langevin 
equation. 

1. Data Collection and Order Parameters Used 


The data are generated from a GSL of the form 
X= a,,Yta,,XYt+ a,3n, + a,4Yn, 


48 


RES mc 


LINEAR ADDITIVE MODEL 


Generated 


Fitted 


40.0 
X(T+AT) 





Figure 6.3 Generated and Fitted Lagrangians for Linear Additive Model. 


Y = a,,X+a,,.XY+ a,,Xn, + a54N, - (6.8) 


Where cross terms are present in the drift (mean) and the diffusion (variance). The 
multiplicative noise is present in the coefficients of the n’s. The APL program 
LANCHESTER described in Appendix D was used to generate the data. Sample 
trajectories of X and Y are shown in Figures 6.6 and 6.7. Figures 6.8 and 6.9 plots the 
sample distributions. 

We now derive the form of the Lagrangian corresponding to the GSL in 
equation 6.8 . | 


The theoretical Lagrangian, L, to be fit to the data is defined in equation 6.9. 


L = 1/2MM- gh ye yMY-g”) , 
ae) 
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Figure 6.4 Generated and Fitted Lagrangians for Non-Linear Additive Model. 
gh = H+ i/ag’, Ogham” , 
gh¥ = gh gy. . (6.9) 


where we have assumed V=0 and ft = a,,Y + a,,XY. ‘Therefore, g! is given by 
equation 6.10. 


dg! Gg! Gg! Gg! 
i= f + 1/29), — ioe ——— 6.10 
g ose a See 2B 1/28" ay (6.10) 
bb <s = , 2 il a 
By = 83,82 = ay’, By = a3%, 8°, = aay, (6.11) 
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Figure 6.5 Generated and Fitted Lagrangians for Linear Multiplicative Model. 


gi = ay,¥ ate ay oXY + 1/2a) yao, (6.12) 
For g* we have 
g? = ee 1/2! an nee oe 1/28, Lat ae ea ’ (6.13) 
and f* = a5)X + a5,XY. Therefore, 
gemma X tea, XY + 1/24, ,a,, (6.14) 


Next we calculate the gi” terms. 


S| 


TRAJECTORY OF X ORDER PARAMETER 


Figure 6.6 Trajectory of Order Parameter X. 


g 8) + B78 = ay + ayy Y? 

BB) + BiyB) = A434) 4X tay yayyY ; 
B18) + 898 78 
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8 > oe 2 
B18", + B82 = a3X” + ayy 


a 





(6.15) 


(6.16) 


(6.17) 


(6.18) 


TRAJECTORY OF Y ORDER PARAMETER 





Figure 6.7 Trajectory of Order Parameter Y. 


Next we need Suv: The g#Y form a 2 by 2 matrix whose inverse is Env: 


Thus 
92? . gl? 
ny = (det eee (6.19) 
. g?! gil 
Ieee? et 2? 


where det ghY = g''g gg * = (det ae . The Lagrangian is 


L=1/2(X-g!)?g,, + 1/2(X-g!)g) (¥-8”) (6.20) 
+1/2(¥-g)g9)(X-g!) + 1/2(Y-B°)’59 


= 1/2(det glY) 1 X-g!)g?!-2(X-g!y ¥-g2)g!2 + (¥-g2)?g! (6.21) 


a3 
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Figure 6.8 Distribution of Order Parameter X. 


since all the terms contain the (det gli factor and since g!? = g?! . 


ive short time conditional probability distribution Ly 
= P(X(t+ At)[ X(t + (n-1)At)) is 
PS = g hr? (2nAt)!/? exp(-L_At) (6.22) 
where 


g. = det (2nwn ; 


= (1/2At)(ME (t+ (n+ LAt)-MM(t+ ndt)-g!At) (6.23) 
X gu y(MY (t+(n+ 1)At)-M%(t+ nAt)-g7At) 


Sd 
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Figure 6.9 Distribution of Order Parameter Y. 


and specifically, 


ee bip22erl2 21) -] ae , J Ze 
Pees 6 eg ge) (CX, “Neg At)’g 
-2g!4x F - X -glAtyYy,” - Yo -g7At+(¥,” - Y,-At)g!!] 
where Xe = X (tt(nt+ I)At), X, = X(t+nAt), and where similar equations exist 
forte Y . 


2. Nlaximum Likelihood Fit of the Lagrangian 
After a functional form of the Lagrangian has been derived (or guessed), the 
next step is to estimate the parameters in the Lagrangian. This is done using a 


maximum likelihood fit, that is, we wish to maximize 
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Tea 
M = Il Ree 


i=l j=1 4 


where Pi = P(X(t+(j+ 1)At)/X(t + jAt)) is the short time conditional probability 
distribution for experiment i, at time t+(j+1)At (post-point). As before, we will take 
the logarithm and pull the minus sign out front, and then minimize the resulting 
function. This becomes 


In M 


» ys In Pi 


p % (i/2in (2nAt)-1/2Ing;;+ L,At} = -N 


We minimize N with respect to the parameters/coefficients of the Lagrangian form. 
This will be a minimization problem of 8 parameters for our example. N will likely 
contain many minima but we are only interested in the best fit, 1e. the global 
minimum. Therefore we will use the simulated annealing algorithm developed for this 
purpose of performing the best fit. Again we are looking for the best fit of the 
Lagrangian to the theoretical Lagrangian and not attempting to extract any particular 
parameters/coefficients. However, any constraints or conditions imposed on the 
parameters/coefficients based on phenomenological reasons should and must be 
included to obtain the best fit. 
3. Results of the Fit 

In addition to the above example, one other fit to data obtained from a 
different GSL is given to test the code and to examine the similarities between the 
Lagrangians. The results are given in Tables 2 and 3. Again we will normalize the 
Lagrangians. For the two dimensional case we will use X; = X_ + ghX (X_) where 
gh are the diagonal terms of the metric and are a measure of the curvature in the 


order parameter space. The renormalized Lagrangian becomes 
The generated Lagrangians are plotted in Figures 6.10 and 6.12. The fitted 


Lagrangians are plotted in Figures 6.11 and 6.13. Again we find good agreement in the 


drift terms or the location of the minimum of the Lagrangian. These are located near 
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TABLE 2 


TWO ORDER PARAMETER MODELS 


Model GSL 
I Linear Additive A=a,Y+ an), + 4)3N, 
BS aap a5 * 255M 
II Non-Linear Multiplicative X=a,,Y+a,,XY+a,,n, +a,,YN, 
pet ioe Mia a4) 
TABLE 3 


TWO ORDER PARAMETER MODEL RESULTS 


Model Generated 
Coefficients 


I an) = -0.008 
a), =0.3 
a,3=9.1 
a5) = -().004 
a,,=0.1 
a,,=0.3 


Il a= 0.01 
a)>= -0.004 
a)3=9.5 
a)4— 9.01 
as) = ().02 
a5, =-0.001 
a5, = 9.015 


a5,=9.7 


7 


Fitted 
Coefficients 


A), = -0.0121839 
4, = 0.0787236 
4), = 0.130564 
A, = -0.00655852 
,, = 0.152455 
4,, = 0.0209012 


A), = 0.0550272 
A, = -0.0073603 
4) = 0.246917 

4 = 0.00235102 
A, = 0.0855577 
a, = -0.00449816 
4, = 0.00857855 
4,4 = 0.281460 


the contour labeled 0.01 in the plots. The contour level spread is a measure of the 
variance and it is evident that there is more spread in the levels in both of the fitted 
Lagrangians. This is somewhat to be expected when dealing with a small sample from 
a probability distribution. More work.is being done here to extend the results to more 
complicated Lagrangians, using more data, and using the simulated annealing program 
to find a better minimum (by including more points in the search). These early results 


have been encouraging. 


_ LINEAR ADDITIVE MODEL:GENERATED 





Figure 6.10 Generated Lagrangian for the Linear Additive Model. 


D. TWO ORDER-PARAMETER MODEL USING JANUS DATA 
1. Selection of Order Parameters 
To ensure a simple description we will assume here that our order parameters 
are the numbers of personnel in each force and that we will look at the attrition, 


similar to a Lanchester approach. 
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Figure 6.11 Fitted Lagrangian for the Linear Additive Model. 


2. Data Collection 

JANUS was chosen as a source of data because of its easy accessibility 
(present at NPS on a number of computers) and the combat data is in an easily 
reducible form for analysis. The high degree of graphics support also gives us a broad 
overview of the scenario we are dealing with. A Blue defense versus a Red Offense 
scenario was used. Terrain was fictitious and generated using the graphics capability of 
JANUS. It was flat with no foliage nor any built-up areas (cities). The defender had 
the capability to perform pop-up maneuvers, 1. e. moving from full defilade to partial 
defilade. This was a separate condition from the terrain. The attacker moved across 
the terrain exposed. Blue was composed of 3 task forces of 12 M6QA3 tanks and 4 
M901 TOW weapons each and was in a line defense. Red was composed of 3 task 
forces of 30 T72 tanks and 9 BMP troop carriers each and was in the attack. Red 


a 


<0 


NON-LINEAR MULTIPLICATIVE MODEL:GENERATED 





Figure 6.12 Generated Lagrangian for the Non-Linear Multiplicative Model. 


objective was to penetrate Blue’s defense, and Blue’s objective was to repel Red’s 
attack (mutually exclusive mussions). 

Twenty runs were collected using the batch mode of JANUS. Data collected 
was driven by the attrition process, and each attrition was classified as an event. At 
each event, clock time of the simulation and attrited side was recorded. Data consisted 
of a collection of clock times and status of forces at clock time. 

3 Development of the Lagrangian 

The first step in the development of the Lagrangian is to examine the sample 
paths of the data. This might suggest an appropriate model, such as Model I or Model 
II to begin with. It is probably best to begin with a simple model and move on to 
more highly non-linear models unles the data is known to be complicated. These 


functional forms can be as complicated as you wish, combining trigonometric, 
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Figure 6.13 Fitted Lagrangian for the Non-Linear Multiplicative Model. 


exponential, or polynonual terms. Hlowever, as mentioned before, a Padé approximant 
can approximate many functional forms and thus a ratio of polynonuals would be your 
best first guess. 
4. Performing the Maximuin Likelihood Fit 

Having selected a functional form we must now. estimate the 
coefficients/parameters in the form. This is done by using a maximum likelihood fit as 
described in the two earlier examples and in Chapter 5. This requires the user, in order 
to use the simulated annealing program, to write a subroutine for his cost function, 1.e. 
the function he wishes to maximize/minimize. This should be written in terms of a 
minimization for the program to work. The user must also decide how many points he 
Wants to use in the search, starting temperatures, any constraints on the 


parameters/coefficients he wants to impose, and a starting point. These are relatively 
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easy to incorporate into the program. One rule of thumb that has been found is to 
start with a small number of points in order to check your subroutine and then 
increase to the amount of points you want. The minimum suggested number of points 
for one dimensional problem is 10000. This usually allows the simulated annealing 
algorithm sufficient time to find a good fit. Of course, the more 
dimensions/coefficients you have, the more points you will need in order to get a good 
fit. The minimum cost and the point associated with that cost is the final output. 
These are your fitted coefficients for your Lagrangian. 

With the Lagrangian, we can now test the fit using similarly obtained data, i.e. 
data from the same scenario. If the fit does not seem satisfactory based on some pre- 
selected conditions, then we must return to the development step and try a different 
functional form. This iterative process continues until we are satisfied we have a good 
fit. 

>. Path Integral Representation 

With a satisfactory Lagrangian we can now examine the long term behavior of 
the svstem using recently developed code for the path integral (not available at the 
time of this thesis). By using the path integral code we can numerically determine the 
probability distribution of our order parameters at any time t. For example, an item of 
interest may be when a particular transition point in a battle might occur. We would 
use the path integral to calculate the probability distribution at each subsequent time 
step and then look for evolving trends of that distribution. For example, imagine we 
had a bistable probability distribution, i.e. one with two most likely states, and wanted 
to know how the system evolved from one minimum state to the other. We would 
then keep stepping through the calculation of the path integral until we found a 
noticeable change in the distribution. 

6. Sensitivity Analysis 

By using sensitivity analysis we can test the robustness of our fit and of our 
functional form. This might proceed as follows. First obtain sufficient data from some 
scenario you wish to study, e.g. our JANUS data. Our particular order parameters are 
the number of each vehicle on each side, so we have four order parameters. We decide 
on a functional form and perform the maximum likelihood fit. After having satisfied 
certain conditions, we are satisfied with our fit. We now return to our scenario and 
make a change in a microscopic variable, for example the pop-up rate of the tanks in 


the defense. Using the same functional form of the Lagrangian, we perform a 
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maximum likelihood fit using the newly obtained data. Remembering that we are 
fitting the Lagrangians and not the individual coefficients/parameters, we are not 
concerned with agreement here. We plot the Lagrangians for both the original data 
and the modified data and determine if there is any significant change. We might 
continue this process for differing values of the pop-up rate or turn to a different 
variable, e. g. the initial force level of each side, max range of the weapons, 
manuevering speed, change in the terrain, weather, etc. If the Lagrangian has not 


changed appreciably we can feel confident that our model is appropriate. 


E. SUMMARY 

In this chapter we presented three examples of the method using the one order- 
parameter and two order-parameter models. In each case we have assumed no explicit 
time dependence in the functional form. This was done to keep the examples simple. 
If the multiplicative noise is small compared to a constant noise term, then the 
Lagrangian is approximately quadratic in the post-point variable, and all distributions 
are Gaussian. However if the multiplicative noise is significant, if the mean is higher 
order than linear, then the long time distribution will be non-Gaussian even though the 
short time distribution is Gaussian. This must be emphasized. Otherwise the model 
would numic a Gaussian distribution and be of limited use. 

We have also presented the use of the simulated annealing algorithm to perform 
the maximum likelihood fit. A relatively new algorithm, simulated annealing seems to 


be particularly useful, especially when little is known about the system. 
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VII. CONCLUSIONS 


A. INTRODUCTION 

We have presented in this thesis an alternative to modeling combat which 
incorporates its severely stochastic nature. We have found it is a promising model for 
different types of combat if certain assumptions are satisfied. A relatively large combat 
system 1s necessary to Satisfy this assumption. This should be of battalion size or larger 
for land battles and a carrier battle group or larger to describe the outer air battle of 
naval engagements. 

A methodology has been developed which will allow the combat analyst to derive 
specific functional forms for use in the decision aid PIACA. The necessary 
mathematical theory has been developed and provides the foundation for the 
methodology. A simulated annealing algorithm to perform the maximum likelihood 
fits was developed. A FORTRAN program was written using the alogrithm and is in 
Appendix A. Three examples using the methodolgy and theory have been given with 
mixed results. These results are discussed next. Following that discussion, an outline 
for the development of the decision aid is given. Finally, the significance of the model 


is discussed. 


B. RESULTS 

The “truly-nonlinear” path integral method has been used sucessfully to describe 
such large scale systems as financial markets, the brain, and in nuclear physics 
[Ref. 8,17,18,19,20,21,22,23,24,25,26,27]. This is the” first Sattempt to Wacenaily 
numerically calculate the nature of large scale combat using these methods 
[Ref. 1,2,12]. In the one order parameter example with constant variance, a 
Lagrangian was fit to the data with excellent results. In this case the path of the 
expected value of the order parameter followed the associated deterministic Lanchester 
equation which has been shown to model combat quite well under certain conditions. 
Other features of the one order parameter model was the incorporation of a drift and 
diffusion term of the order parameters. This is an improvement over the simple 
Lanchester approach. 

The two order parameter example was the first look at incorporating 


multiplicative noise into the model. This is a significant improvement over the 
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Lanchester approach, stochastic or deterministic. Finally, using as our order 
parameters the level of the Blue and Red forces, a two order parameter model was 
developed for a JANUS simulation. 

Given the time constraints of the NPS program, the purpose of this thesis was to 
lay the formulation for future theses to build upon. The Lagrangian we have 
formulated can be further investigated using this approach, e.g., using the simulated 


annealing program to fit a Lagrangian to a specific scenario. 


C. DEVELOPMENT OF THE DECISION AID 

We will now outline below the development of the decision aid PIACA and 
Suggest some design requirements. 

First, to be useful as a decision aid under severe time constraints, PIACA should 
be graphically and qualitatively oriented. The Lagrangian should be presented in a 
form similarly developed in Chapter 4, i.e. a graphical portrayal of most likely states 
and the risks associated with those states. This will allow easy assimilation bv the 
commander and his staff. Of course, human factors should also be incorporated in this 
design. 

Second, PIACA should enable the user to develop his own form of the 
Lagrangian by using data from simulations he has selected. This would allow for 
Lagrangians to be developed which are terrain, scenario, or commander dependent. 

Third, once a Lagrangian is developed, the long time probability distribution, or 
path integral should be calculated easily, 1.e. in real time. This is a heavy requirement 
since at present the path integral is not easy to calculate on supercomputers much less 
something which can be brought to the field. 

To summarize, to have the full decision aid will require: 


e The See of efficient algorithms to solve the foo integral. This will 
most likely be a joint improvement in software and hardware. 


e Graphical devices to portray the evolving long time probability distribution 

e User-friendly software to interface with, the user to Bays/oD) Lagrangians of 
scenarios hé has selected with the capability of pre-selecting Lagrangians from 
Standard scenarios. 

° gem. depiction of alternatives. This is a must requirement for any decision 
aid. 

e Ability to, link to other users. This will then provide an easy information 


exchange in the meta-language of the model, 1.e. passing values of the order 
parameters, scenarios, etc. 
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Obviously, there is much research to be completed in all these areas. One 
purpose of this thesis was to lay the groundwork for the development of the full 
decision aid as well as provide a unified means of describing combat using physical 


concepts. 


D. SIGNIFICANCE OF THE MODEL 

It has been shown that this model is extremely rich in applications. We have 
kept the examples simple to illustrate the use of the method. We have assumed only 
implicit time dependence and a free particle Lagrangian, i.e. V=0. The significance of 
the model could be increased with the addition of these terms. For example, it would 
then be possible to: 

¢ Incorporate boundary conditions, i.e. constraints on the forces either 
geographically or from higher levels of command, bv the addition of the 
potential term, V. Such a model might be useful to describe relatively isolated 
combat. 

e The possibility to examine “phase transitions’, such as during nuclear, 
biological, and chemical exchanges which drastically alter the evolving battle. 
Other relatively small “phase transitions” could be examined such as when a 
commander has reached a critical point in his force level and how he reacts. 
These “phase transitions” could be modeled by matching at the critical point the 
two Lagrangians, i.e. one from the left and one from the right of the critical 
point. 

e The full power of mathematical tools such as stability analysis and calculation 
of first passage times are available for examining the structure of the 
Lagrangian. These tools have been used by physicists and others to examine 
large and complex systems and many results could be applied to combat 
systems. 

Obviously, there is much work to be done in this area that was beyond the scope 
of this thesis. The author hopes more work is done and this thesis layed the 


appropriate groundwork for others. 
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APPENDIX A 
THE SIMULATED ANNEALING ALGORITHM 


Simulated Annealing is a stochastic optimization algorithm for finding global 
extrema. First, the algorithm is described, and a stochastic model, specifically the 
Markov process, is used to show that the algorithm converges to the global optimum. 


Second, the algorithm is used on several unconstrained minimization problems. 


1. BACKGROUND 

We will begin this discussion by first asking the question, what is simulated 
annealing? In order to answer this question, we first need to Know what is meant by 
annealing. Annealing is a process whereby a metal or crystal substance is first melted, 
then subsequently cooled to freezing temperature. At intermediate temperatures the 
substance is allowed to come to equilibrium. The purpose of an annealing experiment 
is to determine the ground or lowest energy state for that particular substance. The 
sequence of temperatures at which the substance is allowed to come to equilibrium is 
refered to as the annealing or “cooling” schedule. If this cooling schedule is too rapid, 
then the state of the substance will most likely be trapped in a metastable state and will 
not reach the ground state and it is said to have been “quenched”. 

Simulated Annealing is then a simulation of this process using a computer. We 
generate an initial configuration of the system and calculate its energy, Ey. Then the 
state of the system is changed, and the new energy, E,, is calculated. If AE = E, - Ey 
< 0, then the state of the svstem becomes this new state. If AE > 0, then the state of 


eAE'T where T is the 


the system becomes the new state with probability = 
temperature. Otherwise, it remains in the old state, and the process is repeated until 
T=0. In this case, we are interested in a cooling schedule which minimizes the total 
energy. The above procedure is refered to as the Metropolis algorithm after 
Metropolis, et. al. [Ref. 28] who developed it to perform equation of state calculations 
for substances at an equilibrium temperature. The algorithm is a useful one for 
performing simulation in statistical mechanics. In large physical systems, one is 
typically interested in the average properties of systems in equilibrium since these are 
the ones that are directly measureable. For example, the procedure has been used to 


simulate ferromagnetism of an Ising model [Ref. 29] and used to calculate <U>, the 
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average energy, and M, the magnetic moment. Also, the entropy of the amorphous 
magnetic state can be calculated using an algorithm similar to the Metropolis 
algorithm [Ref. 30]. 

Kirkpatrick, et. al. [Ref. 31] reformulated the algorithm to be used in the 
minimization of a general class of functions, analogous to the minimization of energy 
in the statistical mechanics system. They used the algorithm to calculate the optimal 
placement of integrated circuit chips on a computer circuit board, and also obtained a 
solution to the N-city traveling salesmen problem. Since that time, the simulated 
annealing algorithm has been used in various forms for problems in circuit design 
[Ref. 31], image reconstruction [Ref. 32,33,34], target tracking, [Ref. 35] layer 
assignment [Ref. 36], speech recognition [Ref. 37], and others [Ref. 38,39,40]. 

In section 2, we discuss the general class of problems associated with non-convex 
optimization. In section 3, the general simulated annealing algorithm is discussed with 
a description of the Markov chain model of simulated annealing. We show that if 
there exists a generation function, and an acceptance function that satisfy certain 
conditions which impose an irreducible, aperiodic Markov chain, then the algorithm 
will converge to the global optimum. In section 4, we give the results of an experiment 
which tested various combinations of generating and acceptance functions on the 
minimuzation of three different cost functions where the global optimum is known. 
Preliminary results of some modifications to the algorithm are also given. We conclude 


our results in section 5. 


2, NON-CONVEX OPTIMIZATION (NCO) 

Non-convex functions are functions which have multiple extrema. In NCO, we 
are typically interested in obtaining the global minimum or maximum. Algorithms 
which attempt to find these extrema can be classified as being either deterministic, 
stochastic or mixed in origin. Some deterministic algorithms such as quasi- Newton 
(BFGS) or steepest descent methods typically locate only local extrema and are only 
guaranteed to find the global when the function to be optimized is convex. There are 
other deterministic methods which have achieved good performance such as Levy and 
Montalvo [Ref: 41] for locating multiple extrema. Stochastic methods include the class 
of algorithms associated with simulated annealing, and which are purely stochastic in 
nature. Mixed algorithms are typically of the multiple starting point, iterative 


improvement variety. An example of this is the IMSL routine ZXMWD which uses a 
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quasi-Newton method with multiple starting points and selects the lowest value. As 
the number of starting points is increased, the probability the global minimum has 
been found is also increased. 

The question now becomes, why use simulated annealing? Obviously, if there are 
algorithms which exploit the structure/properties of a class of problems, then these 
should be used. For example, there exist heuristics for the traveling salesmen problem 
which exploit the characteristics of the problem and probably would perform better 
than using simulated annealing. Simulated annealing becomes very useful if the 
problem vou are interested in lacks anv special structure or for which there are no 
heuristics available. Simulated annealing is easy to implement and gives good insight 
into the problem, if you can identify the cost function to be used. This will become 


more apparent when we describe the algorithm. 


3. GENERIC SIMULATED ANNEALING ALGORITHM 


The algorithm is very simple and is stated as follows [Ref. 42]. 


STEP 0 

e Pick a starting point X,. This could be at random or set to some initial guess. 
If you have some idea of the space to be sampled, the algorithm will run much 
quicker if it is confined to a smaller space. 

e Set the initial temperature, Tj. Again if the sample space is confined, then T, 
can be set at a lower temperature. If not, then set T, to some high value. Xo 
and T, are dependent upon the function to be minimized and are considered 
control parameters. 

e Select a cooling schedule. This will be dependent upon gy discussed in step 1 
below. This is equivalent to setting T(t) = f(t) where t is the step counter. For 


example, set I(t) = T)/(1+t) (inverse linear cooling). 


STEP | 
e Pick a new point, x,. This new point will be picked according to some 
generating function gy which could be a function of the temperature, T. Some 


examples of generating functions are given in section D. 


STEP 2 
e Calculate AC = C, - Cy where C, = f{({x,). This is where the dynamics of the 
cost function enter into the algorithm. 


e Generate a uniform (0,1) random variable, U. 
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e If the acceptance function (discussed below) ay7(AC,T) > U, then accept the 
point x, as the new state of the system and let xX» = x,,Cy = C 

le Cee 
obtained so far) 


1° 


= Cay x = Xp. (This keeps track of the lowest value 


min min 
¢ Otherwise, keep Xp. 


e Repeat step l. 


Another control parameter is the number of iterations/steps to complete. This is 
dependent upon the starting temperature and the number of variables in the function 
to be optimized. Generally speaking, if the sample space is confined to a small region, 
the number of steps needed will be reduced. 

It is clear from the above, how we can model the algorithm as a Markov process 
and use those results to show global convergence. First, we consider the sample space 
as our Markov state space. We can consider this to be a finite state space if we take 
into account the resolution of our computer. Intuitively, it can be seen that if the 
generating function is capable of generating any point in the space, and the acceptance 
function has a finite probability of accepting the point, then every point could be 
visitedgenerated an inifinite number of times. Mathematically, this is equivalent to 
stating that our Markov chain is aperiodic and irreducible. For this to be true, the 
generating function must satisfy four conditions: 

1. gy bea bonafide probability distribution 

2. 87(Xq, X) = Sz_(X), Xy) symmetry. This condition is needed for aperiodicity of 
the Markov chain, i.e. that the generating function is symmetric and no cycles 
are present. 

3. gy(Xq, X,) > O for all X9,x, S (sample space). This could be generated as a 
path from Xp tO X). This says that every state can be reached by any other 
state, although not necessarily in one step. 

4. lim gy(Xp,X,) must exist. 

The three conditions that the acceptance function must satisfy are: 

aT(X9,X, )/az(X Xp) be a monotone positive function. 

2. If C, < C,, then a7(C,,C,) > 0 ie. there must be a positive probability of 

accepting a decrease in the cost function. 


3. lim ay(Xp,X,) must exist. 
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If the aforementioned conditions hold, and if it can be shown that every state can 
be generated an infinite number of times, then the algorithm will converge to the global 
optimum. The preceding results are stated here for completeness of the paper. The 
results are proven in the literature [Ref. 42,43]. 

It can be seen from the above that the simulated annealing algorithm is actually 
a whole class of algorithms depending on vour choice of the generating and acceptance 
functions. For example, for the Metropolis algorithm and many other algorithms in 
the literature, gy is a uniform distribution over a neighborhood of fixed size. I-e. the 
next point generated is some fixed step distance from the previous point akin to a 
Random Walk process. A degenerate case would be to let gy equal some constant, ie. 
the uniform distribution over the whole space. This would be equivalent to a random 
sampling algorithm, keeping the state with the least cost. In our experiments, we use 
generating functions which sample the whole space with a bias towards the current 
point. This is equivalent to a random “hop” process, where the process can “hop” 
around the state space, and hops farther at higher temperatures than at lower 
temperatures. For example, suppose we have some function with multiple hills and 
valleys (extrema). Now suppose we have a ball with a lot of energy (high temperature) 
bouncing about this surface of hills and valleys. The ball loses energy according to 
some friction function, ay, which is dependent on the energy of the ball. As the ball 
loses its energy, it will gradually fall into a valley corresponding to a minima. How the 
ball jumps around and how fast the ball loses energy will determine whether the valley 


reached is the lowest valley. 


4. EXPERIMENTAL RESULTS 

Experiments consisted of testing various combinations of generating and 
accepting functions on three cost functions. The purpose was to determine the most 
effective combination in terms of a measure of performance (MOP). One MOP that 
could be used and has intuitive appeal is CPU time. Another MOP is the acceptance 
ratio which is the ratio of the number of accepted points divided by the total number 
of generated points. These MOP’s are directly related. This can be seen if you 
consider that the times to generate a point for each generating function are essentially 
equal. Since the number of accepted points is equal to the number of steps (step 
counter is incremented when a new point is accepted), and this was held constant, the 


only variable in the ratio is the number of points generated. Therefore, the only 
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difference in the runs was dependent on the number of points generated and so the 
acceptance ratio was used as a MOP. Obviously, this MOP could only be used 
effectively if the algorithm was actually able to locate the global minimum. 


The cost functions used are as follows: 
cl = x4. 16x? + 5x 


This has 2 minima at x=2.90, C=-78.33 (global) and x=2.74, C=-50.06 (local) 
(Ref. &szu2]. 


Cc? = x2. 2y? - .3c0s(37x) - .4cos(4zy) + 0.7 

This has multiple minima with global at x=y=0.0, C=0.0 [Ref. 39]. 
Ci = x2 + 2y? - .3(cos(3mx)cos(4zy)) + 0.3 

This has multiple minima with global at x= y=0.0, C=0.0 [Ref 39]. 
The generating functions used are as follows: 


g7! = (1/m) T(T + x*) Cauchy 


2 2 
gr? = 1//2ne? e% 125° Normal (0, 6? = 1) 


2 2 
1/./2ne7 eX (26 Normal (0, o? 


er 4) 


The Cauchy was chosen because of its infinite variance (wide tails) which should 
indicate a good sampling of the space and not have a tendency to get trapped in a local 
minima early. The normals were chosen to see if the variance had any effect on 
trapping in local minimas. If it did, then that would indicate further research is needed 
to determine if the variance could be used as a control parameter, or if faster cooling 
could be used. 


The acceptance functions used are as follows: 


lI 


a +etAC/T) Heat Bath 
hy* = eAC/T Boltzmann 


The Boltzmann function is used in the Metropolis algorithm and is a 
fundamental function in statistical mechanics. This form arises quite naturally in our 


problem, as our cost function for the parameters is the Lagrangian of the probability 
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distribution of the independent variables. One modification to the above acceptance 
functions for better performance is to link T’, the control parameter, to the relative 
success of the acceptance function, instead of retaining the global temperature, T. That 
is, the temperature that the acceptance function uses would be a different temperature 
than the one tied to the number of steps or generation of points. We will call this the 
2-Temperature generic simulated annealing algorithm. This was developed for this 
thesis in an attempt to obtain better results in higher dimensions of the 
parameter/coefficient space. At this time it is only an ad hoc procedure, but it does 
sausfv the above conditions for the acceptance function. The two temperatures allow 
for a small degree of independence in sampling the space. The global temperature, T, 
controls the sampling of the space by scaling the amount of the jumps in the space. 
The control temperature or acceptance temperature, T’, controls the amount of cost 
difference we are willing to accept at each step, and is subsequently connected to the 
state of the system since the state changes each time we accept a point. The 
FORTRAN program is contained in Figures A.la-A. 1h. 

For each run, Ty = 5, Xp = (10,10) (2Dim), and T(t) = 1/(1 + t) (inverse linear 
cooling) was used. (For 1-dimensional problem, xg = 10). Table 4 gives results listed 


by generating function. Table 5 gives results listed by acceptance function. 


TABLE 4 
RESULTS BY GENERATING FUNCTION 
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TABLE 5 
RESULTS BY ACCEPTANCE FUNCTION 
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First we compare the performance of the generating functions. We notice that 


er did not converge to global optimum in any of the runs, so we discard it. For gr! 


| outperformed gr. Obviously, this is not 


and er? we note that for every run, gy 
conclusive since it is based on a small number of runs. Many different runs with 
different starting points and initial temperatures would be needed for more definitive 


2 outperformed hy! in 


conclusions. As to the acceptance functions, we can see that ar 
every case. Therefore, based on these results, it would seem that the best combination 
is to use the Cauchy as the generating function, and the Boltzmann as the acceptance 


function. 


5. CONCLUSIONS 

This experiment was intended to introduce the reader to simulated annealing and 
to show how it can be used for optimization problems. Further research is needed in 
the area of different generating and acceptance functions, applying the algorithm to 
different cost functions, and extending it to higher dimensions. Preliminary results with 
a normal generating function, with variance as a function of temperature, indicate there 


is More interesting work to be done in this area. 
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6 FORTRAN PROGRAM 

Contained in Figures A.la-A.lh below, is the 2-Temperature version of the 
Simulated Annealing Algorithm FORTRAN program code. The subroutine COSTFN 
is for the 2OP model. User must input the number of dimensions (corresponding to the 
number of parameters/coefficients to be fit and NOT to the number of order 
parameters), number of steps or total number of accepted points he wants generated, 
the starting temperature for the generating function, the starting temperature for the 
acceptance function, the starting point, number of accepted points to print, total 
number of generated points (controling the amount of time the program will run), the 
time increment of the data, the number of runs of the data, and the number of time 
increments. This information must be included as the first line of the data and can be 


unformatted. 


PROGRAM SA 
HII IEE EEE SIE FEE HEHE HE HE IE IE IEE ENE JHE FE FETE FE FE FE IEE IEE IEE JHE IEEE FEF IE IEHE HEHE IEE IEE HEINE IE IE IEE 
*THIS PROGRAM IS A TWO TEMPERATURE VERSION OF SA 
FEI IE IE FE FIE FE IE FE FE IE FE FE IE IE IEE IE HE FEE JE FE FHE IE FEE FEE IE IE JE FEIECIE HE IEIE IE HEHE IE JHE IHE IE IEE HEHE HEHE 
INTEGER NDIM,ACC,NACC ,NTOT,INACC,INTOT,STEPS,TI,T2 
DOUBLE PRECISION T0O,TOA,X0(8),X1(8),XMIN(8), 
1ITEMPG,TEMPA,C1,CO,DELTAC,H,CMIN, 
1LPERACC ,IPERAC ,Z1(1000,21),22(1000,21),DELTAT 
COMMON /COMI/ NDIM,ACC,NACC ,NTOT ,INACC,INTOT,STEPS,TI,T2 
COMMON /COMR/ CO,C1,TEMPG,TEMPA,DELTAC ,»H,PERACC,IPERAC,TOA;, 
1X0 »X1,XMIN,TO,CMIN 
COMMON /LAG/ 21,22,DELTAT »NUMRUN,NTIME 
COMMON /A1/ NMOD 
COMMON /SA1/ MAXTOT 
READ (2,*)JNDIM,STEPS»,T0,TOA>,(XO(N),N=1,8) ,NMOD,MAXTOT ,DELTAT > 
INUMRUN ,»NTIME 
DO 5 I=1,NUMRUN 
5 READ (2,%)} (Z10I,N),N=1,NTIME4+1),(22(015N),N=1,NTIME+41 ) 
c WRITE (35%) (Z1015N),N=1,NTIME+1 ) 
C5 WRITE (35%) (Z22(0I1,N)»N=1,NTIME+1) 
WRITE (3,%) 'NDIM STEPS TO TOA STARTING POINTS NMOD MAXTOT 
1 DELTAT NUMRUN NTIME' 
WRITE(3,* INDIM,STEPS,T0,TOA,(X0(N) ,N=1,8 ) »\NMOD »MAXTOT ,DELTAT » 
INUMRUN ,NTIME 
DO 23 JN= 1,NDIM 
AMIN( JIN) = XO( JN) 
X1(JN) = XO( JUN) 
CALL SIMANN 
STOP 
END 





Figure A.la FORTRAN Program for 2-Temperature Simulated Annealing Algorithm. 


> 


HEFCE IEEE IEE IE HE IE IESE IE JHE DE FEE SE IEDE HE IE IESE IE IE 3EIEIEDE 3E EIDE JE FE IE HEE IE IE 3EIEIE IE IEDE JEIEIEIEIE HEE IE FE FE IEE IEE JE DE 

SUBROUTINE SIMANN 
INTEGER NDIM,ACC,NACC,NTOT , INACC,INTOT,STEPS,TI>,T2 
DOUBLE PRECISION T0,TOA,X0(8),X1(8),XMIN(8), 

ITEMPG ,»TEMPA,C1,CO,DELTAC,H,CMIN, 

LPERACC ,IPERAC,Z1(100,21),22(100,21),DELTAT 
COMMON /COMI/ NDIM,ACC;,NACC,NTOT,INACC ,INTOT,STEPS,TI,T2 
COMMON /COMR/ CO,C1,TEMPG,TEMPA,DELTAC ,H,PERACC,IPERAC;,TOA> 

1X0 »>X1,XMIN,TO,CMIN 
COMMON /SA1/ MAXTOT 
CALL COSTFN(X0,CO) 
CALL INIT 
DO 10 TI = 1, STEPS 
TEMPG = TO/(DFLOAT(TI)) 

20 IF (NTOT.GE .MAXTOT) GOTO 11 

CALL GT 


CALL COSTFN (X1,;C1) 

CALL HT 

CALL PICKPT 

IF (ACC.EQ.1) THEN 
CALL ACCEPT 


ELSE 
GO TO 20 
ENDIF 
10 CONTINUE 
11 PERACC = NACC/REAL(NTOT ) 
WRITE (3,100) NTOT,NACC,PERACC ,CMIN,( XMIN( NN) »>NN=1,NDIM) 
100 FORMAT (1X,'TOTAL GENERATED ',I8,3X%,*'NUMBER ACCEPTED ‘,18,3X, 
1'PERCENTAGE ACCEPTED ‘,F%.2,3X//1X, ‘MIN COST ',E12.6 53x, 
1/1X»>'MIN POINT °,4(E12.651X%)/1X,4(E12.6,1X)) 
WRITE (3,%) ‘TEMPG' ,TEMPG, ‘TEMPA’,TEMPA 
RETURN 
END 





Figure A.lb FORTRAN Program for ee Simulated Annealing Algorithm 
cont). 
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SEIEIEIE IEEE EE HEE HE HEFE HE HE IESE EDEN HE HEE DEE ESE FE HEE EEE FEIESE SESE HEHE FEMHEIEDEDE SE HEE HE FE JE IE IE HE FE HEHEHE DEE HE IE HE 
SUBROUTINE COSTFN(P,COST ) 
DOUBLE PRECISION P(NDIM),COST ,X(1000,21),Y(1000,21) 
1,L1,L2,L3,LD,LX,LY,DELTAX,DELTAY 
INTEGER NDIM,ACC ,NACC ,NTOT,INACC ,INTOT,STEPS,TI,T2 
DOUBLE PRECISION T0,TOA,X0(8),X1(8),XMIN(8), 
LTEMPG, TEMPA,C1,CO,DELTAC,H,CMIN, 
1PERACC , IPERAC,21(100,21),22(100,21) ,DELTAT 
COMMON /COMI/ NDIM,ACC,NACC ,NTOT,INACC ,INTOT ,STEPS,TI,T2 
COMMON /COMR/ C0O,C1,TEMPG,TEMPA,DELTAC,H,PERACC,IPERAC,TOA; 
1X0 »X1 »>XMIN»TO,CMIN 
COMMON /LAG/ X,Y,DELTAT »NUMRUN>NTIME 
COST = 0.0 
DO 10 I 1»NUMRUN 
DO 20 N 1L»NTIME 
DELTAX = X(I,N+1)-X(1I,N) 
DELTAY = Y(I,Nt1)-Y(I,N) 
H1=P( 1 )¥VYCI,N)4PC 2 DeX01T SN YC ISN) +P04 )*P(8)/2.0 
G11=P(3) 
G12=P(4 )*Y(I,N) 
H2=P(5 )¥X(I »N )4+P(6 X01 NYCI N+P 3 )*P( 7972.0 
G21=P( 7 )*X(I,N) 
G22=P(8) 
Q11=G11**2+G12**2 
Q12=G11*G21+G612*G622 
Q22=G21**24+622%*2 
DEN=Q11*Q22-Q12**2 
IF (DEN.LE.0.0) THEN 


COST=1.0E25 
GO TO 99 
ENDIF 


DETG=1.0/DEN 
LX=DELTAX-H1*DELTAT 
LY=DELTAY-H2*DELTAT 
LI=(LX¥¥2 )*Q22/2.0 
L2=-1.0*( LX¥LY )*Q12 
L3=(LY**2 )¥Q11/2.0 
LD=DETG*(L14L2+L3 )/DELTAT 


Cc IF (DEN.LE.0.0) THEN 
Cc ENDIF 
TERM1=-1.0*DLOG( DETG/DELTAT )/2.0 
COST=COST+TERM1+LD 
Cc WRITE (3,%) ‘A',P(1),°B',P(2),'C',P(3),'D* »P(G) 
c WRITE (3,%) "E',P(5),°F',P(6),'G',P(7),'H' ,P(8) 
Cc WRITE (3,%) ‘EXPERIMENT NUMBER',I,'TIME STEP',N 
Cc WRITE (3,%) "DEN OF DETG’ ,DEN, ‘'Q11° ,Q11, ‘'Q12° ,Q12,'Q22° ,Q22 
Cc WRITE (3,%) ‘LX',LX,'H1'>,H1,‘'G11',G11,‘°G12' ,G1l2 
Cc WRITE (3,%) ‘'LY',LY,'H2',H2,‘'G21',G21,‘'G22' ,G22 
Cc WRITE (3,%) ‘'L1°,L1,'L2',L2,°L3',L3 
65 WRITE (3,*) 'LD',LD,'TERM1',TERM1, ‘COST’ ,COST 


20 CONTINUE 
10 CONTINUE 
Cc WRITE (3,%) NTOT,'*GENERATED POINT ',(P(1I),I=1,NDIM),COST 
99 RETURN 
END 


Figure A.lc FORTRAN Program for (eeneetature Simulated Annealing Algorithm 
cont). 
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FEE FEE FEE IEICE FE IE IESE HEDESE SE FE FEE JE SESE SESE IE DE HE FE IE FEE IESE IEEE SESE IESE SE JE IEDE SE SEDE HE SEDME SEE IE FE IESE EEE HE IEE SEE 

SUBROUTINE INIT 
INTEGER NDIM,ACC,NACC ,NTOT,INACC,INTOT,STEPS,TI,T2 
DOUBLE PRECISION T0,TOA,X0(8),X1(8),xMIN(8), 

LTEMPG ,TEMPA,C1,CO,DELTAC »H,CMIN;, 

LPERACC ,IPERAC ,21(100,21),22(100,21),DELTAT 
COMMON /COMI/ NDIM,ACC,NACC ,NTOT,INACC ,INTOT »>STEPS,TI,T2 
COMMON /COMR/ CO,C1,TEMPG,TEMPA,DELTAC,H,PERACC ,IPERAC,TOA;, 

1X0 5X1 ,XMIN,TO,CMIN 
COMMON /7Al1/ NMOD 
INACC 
INTOT 
NACC 0 
NTOT 0 
CMIN co 
TEMPA=TOA 
T2=0 
RETURN 
END 


0 
0 





Figure A.ld FORTRAN Program for (cee Simulated Annealing Algorithm 
cont). 


HIE IE TEI ESE IE IEEE SE HEE DME IE DE FETE JE DE DE SEE FEDE FE IE FEE FETE HEE DE FE DEBE JE DESDE HE DEDE FEBE SEDE FE IEEE DE EE HEE HEHE 
SUBROUTINE GT 
INTEGER NDIM,ACC ,NACC ,NTOT,INACC ,INTOT,STEPS,TI,T2 
DOUBLE PRECISION 1T0,TOA,X0(8),X1(8),XMIN(8), 
LTEMPG, TEMPA 5C1,CO,DELTAC »H,CMIN, 
LPERACC ,IPERAC »Z21(100,21),22(100,21 ),DELTAT 
COMMON /COMI/ NDIM>,ACC,NACC ,NTOT > INACC ,INTOT >STEPS sTI5T2 
COMMON /COMR/ CO,C1,TEMPG,TEMPA,DELTAC ,H »,PERACC ,IPERAC,TOA, 
1XO »X1»XMIN»>TO,CMIN 
DOUBLE PRECISION U,HOP 
REAL WK(3) 
DOUBLE PRECISION DSEED/5959/ 
DO 10 I = 1,NDiM 
CALL GGCAY(DSEED,1,WK,U) 
HOP=TEMPG*U 
X1(I) = XO(I) + HOP 
IF (X1(1I).GT.2.0.OR.X1(1I).LT.-2.0) GO TO 5 
IF (X103).LT.0.0.0R.X1(4).LT.0.0.OR.X1(7).LT.0.0.0OR.X1(8).LT.0.03 
1GO TO 5 
Cc WRITE (3,100) HOP,I,X1(1I),X0(1I),5U 
C100 FORMAT(1X,'AMT OF HOP’ ,5X,E12.655X,'I',2X,I2/1X%>» 'NEWPT ' »5X,E12.6,5 
Cc 15X%»5 ‘OLDPT' »5X,E12.651X%» ‘CAUCHY RANDOM NO';,E12.6) 
10 CONTINUE 
RETURN 
END 





Figure A.le FORTRAN Program for Cae Simulated Annealing Algorithm 
cont). 
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FEE IEEE IEE SESE EE IESE JE IE IEEE SE SE SEE SESE FE IE IE FEE IESE IE IE SE IESE IESE HEE SESE SE SE EE SE SESE 3 IE SESE IE IE IE IE IESE IE HE 3E IE TE 
SUBROUTINE HT 
INTEGER NDIM,ACC ,NACC ,NTOT , INACC,INTOT ,STEPS »TI »T2 
DOUBLE PRECISION T0,TOA,X0(8),XK1(8),XMIN(8), 
1TEMPG, TEMPA,C1,CO,DELTAC,H,CMIN; 
LPERACC , IPERAC ,Z1(100,21),22(100,21),DELTAT 
COMMON /COMI/ NDIM,ACC,NACC,NTOT,INACC,INTOT »>STEPS,TI,T2 
COMMON /COMR/ CO,C1,TEMPG,TEMPA,DELTAC,H,PERACC ,IPERAC,TOA, 
1X0 »>X1,XMIN>,TO,CMIN 
REAL#®8 Y1 
DELTAC = Cl = CO 
Yl = DELTAC/TEMPA 
IF (Y1.GT. 20.0) THEN 
H = 0.0 
GO TO 99 
ENDIF 
IF (Y1.LT.-15.0) THEN 
H = 1.0 
GO TO 99 
ENDIF 
H = DEXP (-Y1) 
C99 WRITE (3,%) ‘DELTAC’ ,DELTAC, ‘TEMPA‘,TEMPA, 'TEMPG',TEMPG, 
C 1‘'Y1',Y1,'H'»,H»,'NTOT' ,NTOT 
99 RETURN 
END 





Figure A.lf FORTRAN Program for ee peramre Simulated Annealing Algorithm 
cont). 


FETE FEI IEF IEE IEE FEE SE SE IE IE IESE HE IE SE IESE IESE FE IESE 3E FEE IESE SESE HE SE IEIE EIEIO IE 3E SESE FE SESE FE IEEE ESE FE DE FEE FETE HE 3E FE 
SUBROUTINE PICKPT 
INTEGER NDIM,ACC,NACC ,NTOT,INACC ,INTOT,STEPS,TI ,T2 
DOUBLE PRECISION T0,TOA,X0(8),X1(8),XMIN(8); 
1LTEMPG ,TEMPA,C1,CO,DELTAC,H,CMIN;, 
LPERACC ,IPERAC ,Z1(100,21),22(100,21),DELTAT 
COMMON /COMI/ NDIM,ACC,NACC,NTOT ,INACC ,INTOT »STEPS,TI,T2 
COMMON /COMR/ CO,C1,TEMPG,TEMPA,DELTAC ,H,PERACC ,IPERAC,TOA, 
1X0 »>X1,XMIN,TO,CMIN 
REAL U 
DOUBLE PRECISION DSEED/6969/7 
CALL GGUBS(DSEED,1,U) 
ACC = 0 
IF (H.GT.U) ACC = 1 
NTOT = NTOT +1 
INTOT = INTOT + 1 
RETURN 
END 





Figure A.lg FORTRAN Program for je oosratre Simulated Annealing Algorithm 
cont). 
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FEI FEHE FEIE FE HE EIE HEE HE EEE HEE IE HE EE EEE HEHE HEHEHE HEHEHE HEHE HIME IIE HE HEHE EHH HIE HEHE HEHEHE HIE HE HE HIHEH 
SUBROUTINE ACCEPT 
INTEGER NDIM,ACC ,NACC »NTOT , INACC ,INTOT ,STEPS,TI,T2 
DOUBLE PRECISION T0O,TOA,X0(8),X1(8),XMIN(8), 
1ITEMPG , TEMPA,C1,CO,DELTAC,H,CMIN, 
1LPERACC , IPERAC »Z1( 100,21 ),Z22( 100,21) ,DELTAT 
COMMON /COMI/ NDIM,ACC ,NACC ,NTOT,INACC,INTOT ,STEPS,TI,T2 
COMMON /COMR/ C0O,C1,TEMPG,TEMPA,DELTAC ,H,PERACC ,IPERAC,TOA, 
1X0 »>X1 ,>XMIN,TO,CMIN 
COMMON /A1/ NMOD 
DO 15 J = 1,NDIM 
XO(J) = X1(J) 
Co = Cl 
NACC = NACC + 1 
INACC = INACC + 1 
IPERAC = REAL(INACC )/REAL(INTOT ) 
IF( IPERAC.GT.O0.5) THEN 
T2 = T2 + 1 
TEMPA = TOA/(REAL(T2)) 
ENDIF 
IF (MOD(NACC,NMOD).EQ.0) THEN 
WRITE (3,%) ‘ACCEPT A POINT' 
WRITE (3,101) INTOT,INACC,IPERAC ,NTOT,NACC, 
ITEMPG ,TEMPA,(X1(N),N=1,NDIM),Cl 
101 FORMAT (1X%,12,1X,I12,1K%,F4.2,1K%,I17,1K%517,3X%>2(E9.3 93% )/ 
11X »4(E9.3,52X%)/1X%,4(E9.3,2X%),E12.6) 
ENDIF 
IF (CO.LT.CMIN) THEN 
CMIN = CO 
DO 17 JJ=1,NDIM 
XMIN( JJ) = X1( JJ) 
ENDIF 
INACC 0 
INTOT 0 
RETURN 
END 





Figure A.lh FORTRAN Program for (RR DEIa Simulated Annealing Algorithm 
cont). 
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APPENDIX B 
THE PATH INTEGRAL 


l INTUITIVE DESCRIPTION 
a. Sum over Paths 
We begin by using our previous definition of the path integral, and 


interpreting the integral as a specific sum over paths. The path integral is defined as 


Ss 


P(X(O)IX(tg)) = fo “JBNexp(- YAtL,) 
n= 
DX = (2mg,7Aty/? TT (2mg_ 2Ary!/2 dx 
n=l 


which is a long time conditional probability distribution of a variable X at some time t, 
given its intial position at tg . How is this a sum over paths? We will look at the one 
dimensional case but the intuitive extension to higher dimensions can be easily made. 

To easily see the path integral description we first look at the random walk 
problem. Suppose we have had a lot of drinks at Tun Tavern (a famous Marine Corps 
establishment of the late 1700’s). Although we do not want to leave, we have an 
inspection tommorrow and need to get home. Tun Tavern is located as X, and home 
is X(t). Being a true drunk, we would have a 50-50 chance of stepping out a certain 
distance in a time increment At. There will be associated with our walk a probability 
of never reaching home (depending on how many drinks we have had!) and thus there 
is a probability distribution of getting home by time t. Suppose we did this every 
night, then each night we would follow a different path home. This is an example of 
Brownian motion. The short time conditional distribution, i.e. the probability of being 
across the street at (z) at time t+ At given you were at Tun Tavern (x) is 


P(across the street,t+ At}|Tun Tavern,0) = 
P(z,ty + Atlx,t,) = (4mDAt)!/? exp(-(x-z)?/4DAt) (B.1) 


where D is the diffusion coefficient or a measure of your drunken state. 


The long time distribution (being at home at time t given you started 
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at Tun Tavern) can be calculated and is 
P(home,t|Tun Tavern,0) = n/./4mD exp[(-x?/4Dt)/</t] . (B.2) 


Suppose again we are fairly drunk but now we have someone pushing us 
according to some rule. Again there is a chance of not making it home, and we want 
to look at our probability of making it home. We will assume our short time 


conditional distribution is given by 
P(X(ty + At|X(ty) = (2mAtg?)!/? exp(-LAt) (B.3) 


where L is the rule used by this person. We again have paths over which we travel 
going from Tun Tavern to home. There are certain probabilities associated with each 
path and we will sum over all the probabilities of the paths to determine our chance of 
making it home. In the limit as our step size becomes continuous and At= 0, we will 
need to integrate over all positions and over all paths and find that our chance of 


making it home is 
P(X(t)|X(tg)= f° * “DX exp (-fLdt) — (B.4) 


which is our path integral. 
b. What does the Path Integral say? 

We can now easily see what the path integral gives us. Ifa particle (or drunk) 
follows a path which is determined by the Lagrangian, we sum over all possible paths 
(which have been weighted by the Lagrangian) and arrive at a probability distribution 
of the particle’s position at some later time. 

The Lagrangian from classical mechanics is T-V where T is the kinetic energy 
and V is the potential. In classical (deterministic) systems there is no path integral 
since there is only one path which is followed with certainty. This is the so-called 


classical trajectory. In classical statistical mechanics the Lagrangian 1s given by 
Li (X= )2/2¢? (B.5) 


where f is the drift and g* the diffusion. 
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2. MATHEMATICAL DESCRIPTION 
a. Quantum Mechanics 
The utility of the path integral is in its ability to arrive at classical mechanics 
as a special case of quantum mechanics in which Planck’s constant fi < < 1 represents 


the noise of the system. In quantum mechanics the path integral is defined as 


K= °° ‘J exp(iS.A)2X (B.6) 


C 
s= f a Ldt (B.7) 


a 


where S is the classical action. 

In contrast to the statistical mechanical classical sum over probability paths, 
in quantum mechanics we compute the sum over the probability amplitudes of the 
paths. There are no paths that are followed with certain probabilities. 

For more information on path integrals in quantum mechanics see Fevnman 
and Hibbs [Ref. 44]. 

b. Statistical Mechanics | 

Statistical mechanics is a branch of physics which attempts to describe the 
relationship between the microscopic properties and the macroscopic behavior. It is 
thus characterized by the investigation of large scale, physical, chemical, and biological 
systems and the search for underlying similarities. A subfield of statistical mechanics is 
concerned with modeling nonlinear systems using the Fokker-Planck equation. This 
equation is the reduced master equation of nonlinear systems and seems to be 
fundamental to many physical systems. Thus its solution would greatly enhance the 
understanding of such systems, including combat systems. 

It is in this context that we have developed the path integral specifically as an 
application of statistical mechanics. For other applications of the path integral in 
Statistical mechanics and applications in general see Schulman [Ref. 9] and Wiegel 
(Ref. 45]. 
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APPENDIX C 
SIMULATIONS/WAR GAMES DESCRIPTION 


1. JANUS 

JANUS is a computer simulation of combat developed by Lawrence Livermore 
National Labratory. It 1s an event-driven, stochastic simulation written in FORTRAN 
and work is underway of an ADA version. It models individual weapons, such as 
tanks, vehicles, helicopters, fixed wing aircraft and personnel as distinguishable entities. 
JANUS currently runs on Digital Equipment Corporation (DEC) mini-computers and 
has Tektromix 4125 color graphics workstations [Ref. 46,47]. JANUS is run from one 
Or more workstations composed of a high resolution graphics terminal, | or 2 graphics 
input tablets with mice, and a DEC VT-100 terminal to communicate with the 
operating system. 

JANUS can be run either in an interactive gaming mode or in batch mode. The 
graphics terminals display the terrain, location of all combat systems under control of 
that workstation and anv enemy units which have been acquired by those combat 
systems. In the interactive gaming mode, a player first plans his operation by placing 
his forces where he wants them to start and can give subsequent movement orders. 
Once the game begins, the player interacts by giving orders to his forces by using the 
mouse and the menu on the graphics terminal. 

In the batch mode, a particular scenario is chosen including an initial plan and 
the computer simulates the combat with no plaver interaction. Results can then be 
captured Via appropriate commands given at the beginning of the run. 

JANUS has a tremendous amount of flexibility in designing any partiucular 
scenario. User has complete control over graphics symbols, weapon system 
parameters, Weapon platform parameters, terrain and weather parameters. 

Currently JANUS is not an optimum simulation to study C? in batch mode 
because of relatively unsophisticated rules for decision-making, e. g., for a line of tanks 
to go around a lead tank stuck in a ditch. 

Currently JANUS 3.2 is installed at TRAC (TRADOC Analysis Center) 
Monterey at their Conflict Simulation Lab (TCSL). They have 4 workstations located 


in 2 adjacent rooms to simulate the Blue and Red forces. 
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2. TWSEAS 

The Tactical Warfare Simulation Evaluation, and Analysis System (TWSEAS) is 
used primarily by the Marine Corps as an instructional tool at the Command and Staff 
College, Quantico, Virginia [Ref. 48]. The computer simulation is designed mainly as 
an aid to assist in the war gaming by providing casualty, intelligence, movement and 
supply reports. It also calculates positions of all forces both enemy and friendly. 
Controllers act as mediator between the computer and the players by inputing tactical 
commands from the players and then reporting to the players results of their 
commands. This can be done over real or simulated communications nets. 

The simulation is stochastic and all results are printed on high speed or PC dot 
matrix printers. This limits the capability of this simulation as an analysis tool. 
Current work is underway to correct this situation and provide an analytic as well as 


training tool. 


3. SOTACA 

The State of the Art Contigency Analysis (SOTACA) is a high resolution 
graphics device combined with powerful decision rule software to provide the 
commander a tool to select contigency alternatives. The software is written in 
FORTRAN and the hardware used is a DEC VAX. It is deterministic and uses as a 
selection process a series of decision rules which are entered by the user. The primary 
means of unit movement are through arcs which are placed by the user. The program 
then selects an optimal routing based on the decision rules and the characteristics of 
the arcs. Although very useful as a planning aid, it was not suitable as a source of 
data because of its non-stochastic nature [Ref. 49]. 

SOTACA was developed for use by the Pentagon and the CINC’s (Commander 
in Chief's, Atlantic and Pacific Forces). 


4. BGTT 

BGITT is the Battle Group Tactical Trainer and is primarily used as a computer 
simulation of the naval warfare environments. Its configuration consists of one control 
and three player stations. At NPS, the stations are subdivided physically in the C? Lab 
by partitions. BGTT runs on a VAX-11/780 with RAMTEK Graphics Display 
Stations and has a high speed printer for paper output of game results. BGTT is used 


primarily as a staff trainer and as an analysis tool for students at NPS [Ref. 50]. 
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5. OTHERS 
a. CARMONETTE 
CARMONETTE is a high resolution computer simulation of combat using 
small unit combined forces. [t is an event driven stochastic simultion which provides 
for intermediate and terminal results, and is used for feasibility studies of alternative 
weapon systems and tactics over varying scenarios [Ref. 51]. 
b. FOURCE 
FOURCE is a deterministic simulation of division level force-on-force combat. 
It is used primarily as a measure of command and control effectiveness in combat 
[Ref. 51]. 
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APPENDIX D 
APL SIMULATION PROGRAMS 


This appendix contains the APL programs LANGEVIN and LANCHESTER. 
They were used to generate data from the Langevin equation and a GSL as described 


below and in the text. 


1. PROGRAM LANGEVIN 
The program LANGEVIN (shown in Figure D.1) generated data from the 


Langevin equation 
X =-1X + en (D.1) 


where n ~ N(0,1) and g=1 (constant variance). The user inputs two vectors, INIT 
and P which are described in the program. The output consists of a matrix where a 
row corresponds to one simulated trajectory of the Langevin equation. The columns 
are the value of the variable (X in this case) at each time increment (t+jAt) where 
j=1,...,T.This data was then used by the simulated annealing algorithm described in 


Appendix A to perform the maximum likelihood fit. 


2, PROGRAM LANCHESTER 
The APL program LANCHESTER (shown in Figure D.2) was used to generated 
data for the 20P example described in chapter 6. Data was generated from the 


following equation: 


X = a,,Y+a,,XY+ a,3n, + a,Yn, 


Y¥ = a)AXta XY+ a.3XN) + ag4No. 


The user inputs, as before, two vectors, INPUT and P, where the vector INPUT 
contains the initial values of X and Y, time increment At, number of experiments, and 
number of time increments to use. Output is in the form of 3 matrices, HISTX, 
HISTY, and EXP. The rows of HISTX and HISTY are the trajectories of the X and Y 
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Figure D.1 LANGEVIN Program. 


order parameters, respectively. The columns correspond to the values of the variables 
at the time, t+jAt, where j=1,. . .,.[T where T=NEXP X DELTAT. The matrix EXP is 
a concatenation of the HISTX and HISTY matrices and was used as input for the 


simulated annealing program of Appendix A. 
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Figure D.2)> LANCHESTER Program. 
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